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In  this  paper,  the  early  stages  of  unsteady  flow  around 
a  cylinder,  started  abruptly  from  rest  and  subsequently  moving 
with  a  constant  velocity  perpendicular  to  its  axis,  are  studied 
by  coupling  discrete  vortex  potential  flow  theory  with  boundary 
layer  theory.  The  flow  is  considered  to  be  laminar  and  incom¬ 
pressible,  and  the  Reynolds  number  is  high.  The  development  in 
time  of  vortex  motion,  characteristics  of  flow  field,  boundary 
layer  separation  point  and  drag  coefficient  are  given  in  detail. 
Results  of  coupled  calculation  contain  the  interactions  between 
boundary  layer,  external  potential  flow  field  and  near  wake. 

The  calculated  development  of  symmetric  vortex  possesses  many 
of  the  observed  features  in  flow  field  visualization.  Pressure 
and  velocity  distributions  are  also  given.  The  drag  coefficient 
obtained  is  in  excellent  agreement  with  experiment.  In  the  calcu¬ 
lation  of  boundary  layer  separation  point,  the  Startford  method 
is  applied  to  the  quasi-steady  boundary  layer  state.  The  calcu¬ 
lated  result  is  better  than  those  given  by  others  and  the  computa¬ 
tion  is  simpler.  Shedding  times  and  locations  of  discrete  vortices, 
which  have  been  generated  before  the  boundary  layer  reaches  the 
quasi-steady  state,  are  determined  according  to  the  unsteady  M.  R.  S. 
criterion.  The  effect  of  these  initial  discrete  vortices  on  the 
main  vortex  development  and  flow  field  is  also  discussed. 


1 


Foreword 


Under  the  laminar  flow  conditions  of  a  high  Reynolds  number,  studies  of 
two-dimensional  unsteady  motion  around  a  cylinder  (started  abruptly  from  rest 
and  subsequently  moving  with  a  constant  velocity)  are  important  theoretically 
and  have  wide  applications.  In  applications,  when  a  missile  moves  at  a  large 
angle  of  attack,  the  characteristics  of  the  discrete  vortex  flow  field  at  various 
lateral  cross-sections  perpendicular  to  the  missile  body  are  quite  similar  to  the 
two-dimensional  unsteady  motion  started  from  rest  around  a  cylinder.  In  computa¬ 
tions  of  nonlinear  aerodynamic  force  with  a  large  angle  of  attack,  these  two 
above-mentioned  motions  are  alike^.  The  two-dimensional  unsteady  motion  around 
a  cylinder  is  the  basis  for  analyzing  discrete  flow  (at  the  leeward  surface) 
characteristics  of  a  missile  body  with  a  large  angle  of  attack.  Studies  of 
vibrations  of  many  tall,  large  industrial  buildings  and  water  conservation 
facilities  (such  as  television  towers,  tall  chimneys,  and  large  bridges)  because 
of  flow,  and  studies  of  atmospheric  circulation  related  to  environmental  contamina¬ 
tion  are  also  closely  related  to  flow  in  this  model.  In  the  theoretical  sense, 
in  the  case  of  high  Reynolds  numbers,  flows  in  an  entire  unsteady  process  are 

i 

very  complex,  including  separation  of  unsteady  boundary  layer,  formation  and 
development  of  vortex,  secondary  vortex,  secondary  vortex  pair,  unsymmetrical 
motion  of  vortex,  and  formation  and  motion  of  Kantian  vortex  streets.  The  theoret¬ 
ical  analyses  should  simultaneously  consider  interactions  of  boundary  layer, 
external  flow,  and  flow  near  wake.  Mathematically,  this  is  a  high-order  nonlinear 
problem.  The  various  theoretical  methods  at  present  are  incapable  of  serving  in 
a  detailed  description  of  the  entire  complex  process.  Boundary  layer  theory 
cannot  be  used  to  study  flow  in  discrete  regions;  the  existing  boundary- layer 
solution  (2-6]  do  not  allow  for  the  mutual  influence  of  boundary  layer  and  external 


1  The  paper  sums  up  two  other  papers,  "The  Divergent  Vortex  Model  of  Unsteady 
Motion  Around  a  Cylinder  and  Study  of  Its  Initial-stage  Motion,"  and  "The  Initial 
Flow  of  Laminar  Flow  Unsteady  Motion  With  High  Reynolds  Number  Around  a  Cylinder." 
The  first  paper  was  read  in  May  1980  at  the  founding  conference  of  the  China 
Aeronautics  Association  Annual  Science  Conference.  This  persent  paper  is  one 
of  First  Asia  Fluid  Dynamics  Conference  Papers,  delivered  in  December  1980. 
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flow  in  an  unsteady  separation  process.  The  convergent  matching  solution  [6-10] 
related  to  the  two-dimensional  N-S  equation  is  convergent  only  for  a  relatively 
short  motion  time  (V.// R  at  high  Reynolds  numbers;  therefore,  the  solution 

can  only  reveal  the  initial  flow  variation  after  abrupt  start. 

An  entirely  numerical  solution  of  N-S  equation  can  give  a  detailed  structure 
of  the  entire  flow  field,  providing  various  approximate  solutions  with  a  standard 
for  comparison.  However,  in  the  case  of  high  Reynolds  numbers,  many  difficulties 
exist  in  the  unsteady-motion  numerical  solutions;  moreover,  the  numerical  solutions 
fail  to  clearly  reveal  flow  variability;  it  is  not  convenient  to  analyze  flow 
variations  caused  by  different  parameter  variations;  Thus,  there  have  been 
scant  results  on  this  aspect,  mainly  described  in  the  studies  [11-16],  Within 
the  range  of  higher  Reynolds  numbers,  the  main  flow  phenomenon  in  the  separation 
zone  near  the  wake  is  the  vortex  motion  and  vortex-surface  variation.  The  experi¬ 
mental  studies  [17]  also  reveal  that  after  an  abrupt  start  and  moving  around  a 
cylinder,  viscosity  has  very  little  effect  on  motion  and  development  of  a  vortex 
in  the  intermediate  stage  of  the  motion;  the  variation  with  time  in  the  values 
of  the  circulation-zone  characteristics  is  not  related  to  the  Reynolds  numbers. 

So  nonviscous  flow  theory  can  be  used  (with  assistance  from  the  concentrated 
vortex  amount)  to  approximately  represent  actual  flow  of  low-viscosity  fluid. 

This  is  the  potential  flow  model.  The  early-stage  concentrated-vortex  model 
[18-20]  can  be  used  to  estimate  drag  variation  of  the  initial-stage  motion  around 
the  cylinder;  however,  after  a  longer  time  the  drag  drops  rapidly,  even  appearing 
as  a  negative  value,  with  irrational  results  of  feedback  by  the  concentrated 
vortex  toward  the  vortex  surface,  pressure  interruption  on  the  object  surface, 
and  a  force  couple  acting  on  the  system  composed  of  vortex  surface  and  vortex 
pair.  The  report  [21]  combines  the  concentrated  vortex  model  and  the  boundary 
layer  theory,  considering  the  interaction  of  unsteady  boundary  layer  and  external 
flow;  variations  of  the  separation  points  of  the  initial  motion  are  obtained  but 
the  calculated  drag  considerably  deviates  from  the  actual  result.  The  concentrated 
vortex  model  does  not  take  into  account  the  formation  and  variation  of  the  vortex 
surface,  as  well  as  its  influence  on  the  flow.  The  model  cannot  correctly  reveal 
the  distribution  and  transmission  of  the  vortex  amount  and  the  complex  variability 
of  flow  in  the  separation  zone.  The  convergent  multivortex  model  improves  on 


the  concentrated-vortex  model  by  approximately  indicating  places  of  concentrated 
or  continuous  distributions  of  the  flow  field  by  many  scattered  point  vortexes 
or  linear  vortexes.  The  discrete-vortex  model  can  relatively  clearly  describe 
the  flow  corresponding  to  a  longer  time  of  motion;  the  required  computations  are 
much  lighter  than  the  numerical  solutions  of  the  N-S  equation.  This  is  a  very 
advantageous  feature  in  studying  unsteady  flow  problems  at  high  Reynolds  numbers; 
these  studies  are  mainly  represented  by  the  reports  [22-31]  and  recently  by  [31] 
and  [28].  These  investigations  have  made  certain  gains  in  describing  the  overall 
results  of  long-time  motion,  such  as  the  periodic  variation  of  lift  and  drag,  and 
the  general  forms  of  the  Karmen  vortex  street.  However,  these  past  studies  did 
not  consider  the  time  of  occurrence  of  the  initial  discrete  vortex  (newborn  vortex 
and  subsequent  discrete  vortex)  and  the  effect  of  position  on  vortex  development 
and  flow.  Methods  of  determining  newborn  vortexes  are  random.  In  addition,  many 
past  studies  were  numerical  studies,  mostly  with  random  assumptions  and  few 
analyses  as  to  the  physical  regime. 

In  a  situation  of  laminar  flow  with  high  Reynols  numbers,  the  paper  uses  the 
method  of  combining  the  discrete  vortex  model  and  boundary- layer  theory  in  studying 
the  regime  of  initial -stage  flow  of  unsteady  motion  around  a  cylinder,  to  describe 
the  variation  with  time  of  boundary- layer  separation  points,  formation  and  develop¬ 
ment  of  symmetrical  vortexes  and  vortex  surface,  distribution  of  pressure  and 
velocity  at  the  object  surface  at  different  times,  as  well  as  variations  of  drag 
due  to  the  cylinder.  Solutions  in  the  paper  consider  the  interaction  of  boundary 
layer,  external  flow  and  flow  near  the  wake  of  the  separation  flow.  The  paper 
uses  the  M.  R.  S.  unsteady  separation  criterion  to  determine  the  detachment 
position  and  time  of  newborn  vortex  and  subsequent  discrete  vortex,  and  to  discuss 
the  influence  on  vortex  motion  and  flow  of  the  discrete  vortex  produced  before 
quasi-steady  variation  of  the  boundary  layer. 

Model  and  Method 

The  paper  studies  high-Reynolds-number,  incompressible  laminar  flow.  The 
main  influence  of  viscosity  is  limited  to  a  very  thin  boundary  layer;  after 
separation,  the  shearing  layer  is  also  very  thin,  as  a  thin  vortex  layer.  The 
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vortex  quantity  comes  from  the  continuous  separation  of  the  boundary  layer  and 
is  fed  by  a  small  discrete  nonviscous  vortex  with  certain  intensity  and  detaching 
continuously.  Discrete  vortices  compose  a  vortex  surface.  The  trend  of  motion 
of  the  discrete  vortex  is  determined  by  the  local  potential  flow  field;  its  locus 
varies  with  time  and  constitutes  vortex-surface  motion.  With  the  lengthening  of 
time  and  continuous  motion  of  the  discrete  vortex  in  the  flow  field,  the  number 
of  vortexes  increases  continuously,  thereby  causing  unsteady  variation  in  the 
flow.  From  the  solution  of  potential  flow,  the  distribution  of  velocity  and 
pressure  at  object  surface  are  obtained  to  then  calculate  the  separation  of  the 
boundary  layer  and  to  determine  the  position  and  intensity  of  newborn  vortexes. 

The  addition  of  newborn  vortexes  into  the  original  vortex  group  forms  a  potential 
flow  field  of  the  next  interval;  this  result  can  be  used  to  calculate  new  separation 
points  and  correspondingly  to  determine  new  discrete  vortexes.  By  these  repeti¬ 
tions,  the  entire  flow  field  is  composed  of  a  nonseparated  flow  around  a  cylinder 
and  discrete  vortexes  and  their  mirror-image  vortexes.  The  above-mentioned  model 
considers  the  interaction  of  the  boundary  layer,  external  flow  and  flow  near 
wake  in  the  separation  flow  process. 

Nonviscous  Flow  Field 

The  unsteady  motion  around  a  cylinder  (started  abruptly  from  rest  and  subse¬ 
quently  moving  with  a  constant  velocity)  is  equivalent  to  the  abrupt  addition  of 
an  unsteady  motion  (around  a  cylinder)  to  a  uniform  flow.  This  paper  neglects  the 
effect  on  fluid  by  abrupt  placement  of  a  cylinder.  Except  for  the  portions  of 
boundary  layer  and  the  vortex  nucleus,  the  circulation  flow  field  can  be  considered 
as  a  nonviscous  flow  field,  whose  velocity  potential  satisfies  the  Laplace 
equation. 

Definitions:  dimensionless  composite  potential  w  —  w/RV„,  velocity 
potential  <j»  —  $/RVm,  flow  function  <t>  —  4>/RVm,  circulation  quantity 
r  —  t/ZxRVm,  velocity  V » time  /— P.J/R..,  pressure  coefficient 

e,—(p  —  pm)  /—  pV1.,  drag  coefficient  cD  —  D/  ~^pVi.R,  and  lift  coefficient 

eL  —  Li  i  pV'mR.  5 


(-v  indicates  the  corresponding  dimension  quantity,  and  R  and  V  represent, 

00 

respectively,  cylinder  radius  and  velocity  of  incoming  flow). 


The  forms  and  the  flow  equation  and  boundary  condition  in  polar  coordinates 
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The  composite  function  composed  of  the  incoming  uniform  flow  passing  around 
a  cylinder  and  N  vortexes  outside  of  the  cylinder  can  be  derived  by  the  circle 
theorem  [32].  In  polar  coordinates,  the  expression  equation  of  the  composite 
function  is 

H^C*)  —  *  +  —  —  *r.k>  (*  —  «.)  +  y)  «T.tn  (* - 

*  «-i  .-t  v  -*./ 

y  .v 

~  X!  »T,lna  4-  X,  iT.lnQ'Q  (2) 

•» 1  •» 1 


In  the  equation,  zn  is  the  position  of  the  n-th  point  vortex  as  a  function  of  t; 

Tn  is  the  intensity  of  the  n-th  discrete  vortex;  Tn  does  not  vary  with  time 
after  its  determination,  and  the  positive  sign  indicates  a  counterclockwise  vortex. 
During  the  detachment  time  of  separation  of  each  discrete  vortex,  the  potential 
flow  field  should  satisfy  Kelvin's  theorem  of  conservation  of  vortex  quantity. 
Hence,  there  should  be  superposed  on  a  vortex  quantity  with  an  intensity  the  same 
(but  opposite  in  direction)  as  the  detaching  discrete  vortex  on  the  cylindrical 
surface.  After  neglecting  the  integer  items  not  related  to  z,  the  final  simplifi¬ 
cation  form  of  the  composite  function  is 
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The  imaginary  and  real  parts  of  the  above  equation  can  be  used  to  derive 
the  potential  function  and  flow  function  as  the  following: 
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The  velocity  of  flow  field  under  the  polar  coordinates  is 

,  a.  r1  —  1  a.-Vr  £a-£°  ~ 

K,(r,0)  —  — —  co*0  +  2_i  r* - 5 


R. 


_  V  r  ^  sa  ~zJ!.al 


(6) 


K9(r,9) 


-'i±-1  *.«-  ±r.rr--r'c°<9--J 

r-  .-I  K» 

.  V1  r  f  ~~  r-co»  ~ 

.Tt  *  Rj 


(7> 


The  velocity  distribution  at  the  object  surface  is 
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In  the  equation. 
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rr,  sin  9  —  sin  9. 
rr«cos0  —  cos0« 


i  ^  sin  0  —  r »  sin  8, 
rcosd  —  r.cosfl. 


1  +  —  2rr„co*(0  —  0.) 

r‘  +  r',  —  2rr,cos(0  —  0.)  1 
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The  velocity  of  point-vortex  motion  is  the  velocity  of  the  local  flow  field; 
the  motion  locus  can  be  derived  from  a  step-by-step  integration  of  velocity. 

During  the  computation,  the  influence  of  the  induced  velocity  by  the  vortex  itself 
should  not  be  included.  The  vortex  of  the  k-th  point  vortex  in  polar  coordinates 
should  satisfy  the  following  equations: 
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In  the  equation,  S\—R, (rk,8i,r.,t),'),  s7  —  R,(.rk,dk,r.,d,),  ■> 

(n  -  1,2,---„V:  *-l,2,---„V)  I  <12> 

F  should  be  determined  by  the  vortex  quantity  at  the  separation  point  of  the 
boundary  layer.  With  the  2N  differential  equations  mentioned  above,  the  motion 
state  of  N  point  vortices  can  be  determined. 


By  Lagrange  integration,  we  find  that  the  pressure  distribution  at  object 
surface  of  unsteady  motion  (neglecting  the  mass  force)  is 

c,(0,*)  -  I  -  v\  +  2  2  r.  £  «?*'  -  ( 13  > 

Noticing  that  jw  —  dW,Q&-  —  —  — .  on  the  cylindrical  surface,  while 

9»  di  dt 

computing  the  item-by-item  residual  value  to  and  ~q~ 

finally  the  pressure-difference  drag  and  lift  of  the  cylinder  are  derived: 

c°~  ~4*S r*  i(1  +  7i)^sine*  +  (r»"7;N)^co,94  (14> 

<tW 

In  the  equations,  the  definition  of  the  composite  velocity  is  ^  “  «  — >v. 
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Boundary  Layer  Separation  and  Intensity  of  Discrete  Vortex 

At  high  Reynolds  numbers,  the  flow  closely  in  contact  with  the  surface  of 
an  object  should  satisfy  the  unsteady  two-dimensional  incompressible  boundary 
layer  equation  and  the  corresponding  initial-edge-value  condition.  The  outer- 
edge  condition  of  the  boundary  layer  is  the  potential  flow  velocity  at  the  object 
surface  mentioned  above.  The  solution  to  the  internal  and  external  layer  flows 
should  be  obtained  by  coupling;  therefore,  the  problem  can  be  reduced  to  solving 
for  the  unsteady  boundary- layer  separation  points. 

Previously,  all  studies  of  unsteady  motion  around  a  cylinder  by  the  discrete 
vortex  model  used  the  Prandtl  separation  criterion  to  determine  the  separation 
points  of  the  boundary  layer.  At  t=0.35  after  starting  flow  motion  around  a 
cylinder,  separation  begins.  First,  separation  points  appear  at  the  rear  station¬ 
ary  point;  however,  for  unsteady  motion,  separation  does  not  certainly  occur  at 
T=0  [33].  Moore-Rott-Sears  proposed  the  M.R.S.  criterion  [34]  for  determining 
the  boundary  separation  of  unsteady  flow.  The  clear  meaning  of  "separation" 
indicates  the  distortion  of  boundary  layer  characteristics;  the  boundary  layer 
becomes  very  thick  with  a  singularity  appearing  in  the  solution,  whether  or  not 
it  is  the  beginning  of  separation  bubbles  or  the  wake.  According  to  this  crite¬ 
rion,  generally  the  separation  points  are  not  on  the  object  surface.  When  the 
unsteady  motion  tends  to  exhibit  constant  variation,  the  M.R.S.  criterion  is 
simplified  into  the  Prandtl  criterion.  After  an  abrupt  start  of  an  unsteady  motion 
around  a  cylinder,  the  first  separation  occurs  at  t=0.65  and  the  position  is  in 
the  neighborhood  of  40°  from  the  rear  stationary  point.  There  are  appreciable 
differences  of  the  positions  of  separation  points  at  the  beginning  of  the  motion 
during  computations  of  the  two  separation  criteria;  after  t?l,  these  two  criterion 
computations  gradually  approach  each  other.  After  a  further  longer  time,  there 
are  no  practical  differences  [35]  between  these  two  computations.  If  the  separa¬ 
tion  and  motion  of  discrete  vortex  can  approximately  represent  vortex -quantity 
transmission  after  the  separation  of  boundary  layer,  especially  in  studies  of 
flow  after  abrupt  start,  separation  points  derived  by  the  M.R.S.  criterion  should 
be  relied  on  in  determining  the  detachment  time  and  position  of  the  newborn  vortex 
and  the  subsequent  discrete  vortex,  so  the  discrete  vortex  model  can  be  established 
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on  a  foundation  with  clearer  physical  meaning.  Because  the  determination  of  this 
method  of  separation  points  involves  a  numerical  method  to  determine  the  process 
of  boundary- layer  singularity,  this  paper  does  not  repeat  this  complex  computation. 
In  considering  the  influence  on  flow  by  a  discrete  vortex  produced  not  long  (tsl) 
after  the  abrupt  start  of  the  motion,  the  paper  directly  quotes  the  result  given 
by  the  reports  [35]. 

In  investigating  the  variation  of  boundary -layer  separation  points,  notwith¬ 
standing  the  incessant  variation  of  the  external  flow  field  with  time,  the 
variation  of  separation  points  only  requires  a  short  time  of  approximately  t=l 
to  become  quasi-steady  variation  [28  and  31].  In  this  case,  the  corresponding 
separation  points  have  moved  by  approximately  70°  to  positions  in  the  neighborhood 
of  the  rear  stationary  point.  Later,  the  boundary -layer  separation  can  be  computed 
in  a  quasi-steady  state.  This  paper  applies  the  Startford  two-layer  model  method 
[37]  for  computing  boundary- layer  separation  points  to  the  quasi-steady  boundary 
layer.  Noting  that  the  various  physical  quantities  are  a  function  of  time  and 
simultaneously  considering  that  the  extreme  pressure  and  velocity  values  at  the 
object  surface  during  unsteady  motion  are  possibly  not  at  the  same  position,  at 
the  separation  points,  the  following  relationship  should  be  satisfied: 

[4 *'£)']  -001”'  I 


In  the  equation,  p  .  is  the  minimum  pressure,  V*  is  the  object-surface  external 
fringe  velocity  corresponding  to  the  minimum-pressure  point,  and  0BL  is  the 
starting  point  of  the  Blasius  flat  plate  boundary  layer  used  in  the  boundary- 
layer  two-layer  solution  in  the  Stratford  method,  the  conditions  of  two-equal- 
momentum  thickness  at  the  minimum  pressure  point,  the  boundary  layer  of  actual 
flow  and  the  imaginary  flat  plate  boundary  layer,  are  used  to  determine  the 
starting  point.  The  so-called  coupled  solution  is  the  simultaneous  solution  of 

equation  (16)  and  the  external  flow  field. 
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The  intensity  of  the  discrete  vortex  is  equal  to  the  outgoing  vortex  flow 

of  the  boundary  layer  at  the  separation  points  in  a  time  step  interval.  It  is 

assumed  that  in  a  unit  time,  the  vortex  flow  at  the  boundary- layer  separation 

points  is  tT  ;  the  circumferential-direction  velocity  in  the  boundary  layer  is 
3/ 

Vq,  and  the  potential  flow  velocity  at  the  separation  point  is  Vgg,  then 


dr 

’  f*  dt>.  ,  1 

dt 

_ 

~  2* 

11  o7Vadr\. 

(17> 

r. 

_  l_ 

4* 

V\^ 

(13> 

In  the  equations.  At  is  the  time  step  interval.  The  above  equations  indicate  that 
no  direct  connection  exists  between  the  discrete-vortex  intensity  and  the  flow 
details  in  the  boundary  layer;  the  discrete-vortex  intensity  is  only  related  to 
the  potential  flow  velocity  at  the  external  edge  of  the  separation  point. 

Brief  Description  of  Computations  and  Main  Results 

The  work  in  the  paper  includes  the  following  aspects:  the  first  is  the 
application  of  discrete-vortex  model  to  investigate  the  flow  structure  of  unsteady 
initial-stage  motion,  i.e.,  the  occurrence  and  development  of  the  vortex  and 
the  variation  of  flow  characteristics  with  time.  The  variation  rule  of  the 
boundary  layer  separation  point  can  refer  to  the  reports  [35],  [28]  and  [11]  as 
expressed  with  an  empirical  relationship  equation,  such  as 

0,  -  98.0  -  58.0exp[  —  1. 1968(/  -  0.65)]  (19)' 

When  t«l  in  the  above  equation,  corresponding  to  separation  points  determined 
by  using  the  M.R.S.  criterion,  computational  values  are  derived  after  a  longer 
time  approaching  the  criterion  of  stress  at  the  zero  wall.  The  above-mentioned 
computations  are  briefly  called  Section  A  computations,  and  A  is  marked  on  the 
curve.  The  second  aspect  is  the  coupling  solution  of  the  initial -stage  separa¬ 
tion  flow  field  by  using  the  models  and  methods  described  in  the  paper.  This 
part  of  the  computations  includes  two  situations:  the  first  situation  begins  to 
be  computed  when  the  boundary  layer  variation  tends  to  a  quasi -steady-state 
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situation  and  the  positions  of  the  newborn  vortexes  are  given  by  the  equation 
(19)  as  they  are  situated  at  a  place  +70°  from  the  rear  stationary  point  and 
the  detachment  time  is  tQ=1.2.  By  using  equations  (8) -(18)  to  make  coupling 
computations,  the  second  situation  takes  into  account  the  effect  on  flow  of  the 
discrete  vortexes  produced  during  separation  at  t<1.2.  The  time  when  computa¬ 
tion  begins  is  tQ=0.65  and  the  position  of  the  newborn  vortex  is  situated  at 
0°=+4O®.  After  t=1.2,  the  flow  field  coupling  computation  is  carried  according 
to  equations  (8)  to  (18).  Computations  of  the  two  corresponding  situations 
mentioned  above  are  called  computations  of  sections  B  and  C,  and  curves  in  the 
diagram  are  represented,  respectively,  by  B  and  C.  This  paper  uses  the  Runge- 
Kutta  method  to  compute  the  vortex  motion  locus,  and  the  Simpson  integration 
and  proportional  root  solution  method  to  compute  the  separation  point.  For  each 
time  step,  two  discrete  vortexes  are  added  and  at  the  same  time  six  equations  are 
added,  not  including  the  number  of  equations  serving  in  computing  the  separation 
points.  In  the  entire  process,  the  length  of  the  computation  step  is  0.125;  for 
comparison,  the  step  length  of  partial  computation  is  0.1  and  0.2.  For  computations 
of  section  A,  when  t=3.9  each  vortex  surface  is  composed  of  27  discrete  vortexes; 
totally,  162  equations  should  be  simultaneously  solved.  The  initial  position 
of  the  discrete  vortexes  should  have  a  distance  from  the  cylindrical  surface  not 
greater  than  the  thickness  of  the  local  boundary  layer.  For  computation  conve¬ 
niences,  refer  to  thickness  computation  of  cylindrical  boundary  layer,  by  taking 
r=1.05.  The  closing  of  the  newborn  vortex  from  the  separation  point  will  move 
forward  faster  the  positions  of  separation  points.  In  order  to  reduce  this  effect, 
the  paper  adopts  a  method  of  delaying  a  step,  i.e. ,  after  producing  a  newborn 
vortex  for  an  interval  of  a  time  step,  its  inducing  function  begins  to  take  into 
account . 

The  development  variation  of  vortex:  Fig.  1(a)  to  (d)  describe  the  computa¬ 
tion  results  of  section  A  of  vortex  motion.  With  the  lengthening  of  time  and 
the  development  of  separation  at  the  corresponding  boundary  layer.  Fig.  1  gives 
the  formation  and  roll-up  process  of  vortex  surface  behind  the  cylinder,  giving 
the  development  of  symmetric  vortex  of  spiral  shape.  The  motion  trend  of  discrete 
vortexes  is  rational  and  the  vortex  surface  is  of  sequential  shape.  The  vortex 
surface  is  naturally  composed  of  moving  discrete  vortexes,  not  requiring  the 
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method  of  re-separation.  The  computation  result  in  the  paper  is  very  similar 
to  the  figures  and  images  shown  by  conventional  experiments.  In  (a)  through  (d) 
of  Fig.  2,  the  computation  results  of  section  B  are  shown.  Compared  with  the 
results  of  section  A,  the  positions  of  vortexes  in  section  B  are  higher,  position 
ing  at  the  shoulder  part  behind  the  cylinder  with  narrower  flow  region,  flatter 
vortex  shape,  and  vortex  roll  not  as  sufficiently  developed  like  that  in  section 
A.  Hence,  the  different  detachment  time  and  positions  of  discrete  vortexes  at 
the  initial  stage  will  influence  the  variation  of  vortexes  behind  the  cylinder. 
There  is  some  randomness  in  handling  the  separation  time  and  position  of  the 
initial  vortex  in  the  reports  [21],  [24],  [SO]  and  [51].  The  study  [30]  consid¬ 
ers  that  the  position  of  the  newborn  vortex  is  not  important;  this  view  does 
not  have  a  solid  basis.  As  revealed  by  computation,  if  the  shedding  positions 
of  newborn  vortexes  are  not  rational,  or  interruptions  occur  in  the  subsequent 
step  length  while  detaching  of  discrete  vortexes,  or  the  velocity  of  forward 
motion  of  separation  points  varies  greatly,  this  will  induce  irregular  variation 
of  shapes  of  vortex  surface. 

Figures  3  and  4  show  the  distribution  of  pressure  and  velocity  at  the  object 
surface  (with  time)  in  the  flowing  process.  From  the  figures,  the  flow  character 
istics  at  the  cylindrical  surface  when  tsl  are  little  different  from  the  distribu 
tion  of  potential  flow  without  separation.  With  increase  in  time,  because  of 
lengthening  and  effect  of  vortexes  behind  the  cylinder,  when  t=1.4  two  new 
stationary  points  begin  to  appear  in  the  rear  half  part  of  the  cylindrical 
surface,  i.e.,  the  appearance  of  a  very  small  circulation  flow  zone  and  gradually 
subsequently  expanding;  one  terminal  includes  the  rear  stationary  point  and 
another  terminal  extends  forward,  gradually  becoming  a  stabilized  position.  In 
the  neighborhood  of  the  stationary  point,  there  is  a  low-pressure  zone;  its 
range  is  expanding,  corresponding  to  the  increase  in  time.  The  value  of  the 
low-pressure  peak  may  exceed  the  minimum  pressure  value  at  the  front  half  of 
the  cylinder,  and  the  maximum  velocity  at  the  object  surface  in  the  circulation 
flow  region  also  may  exceed  the  velocity  of  the  incoming  flow.  The  peak  value 
of  pressure  and  the  position  of  the  extreme  velocity  value  do  not  coincide 
because  of  the  unsteady  effect.  Figure  3  shows  a  computational  example  of  a 
coupling  computation  of  pressure  distribution;  from  the  figure,  computations 
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Fig.  1  (a)  to  (d).  The  formation  and  development  of  vortex  (A). 

at  sections  (B)  and  (C)  show  greater  differences  only  in  the  neighborhood  of 
the  rear  stationary  point.  This  is  because  curve  C  takes  into  account  the  effect 
of  discrete  vortexes  produced  prior  to  t*1.2;  in  comparing  curves  C  and  A,  we 
see  that  these  two  are  relatively  close.  The  computations  reveal  that  when  the 
discrete-vortex  motion  is  irregular,  and  the  point  vortexes  are  very  close  to 
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Fig,  4,  Velocity  distribution  at  object  surface  (A) 


Fig.  6.  Variation  of  minimum  pressure  points. 
Key:  (a)  Numerical  solution  [11]  of  N-S  equa 
tion;  (b)  Results  of  this  paper  A. 


the  object  surface,  computation  of  pressure  or  velocity  departs  from  the  actual 
cases  because  of  the  singularity  effect  of  the  induced  velocity.  After  passing 
through  peak  value,  the  circulation  flow  along  the  object  surface  undergoes 
velocity  reduction  and  pressure  rise  within  a  considerable  period  of  time; 
attention  should  be  given  to  the  new  flow  phenomenon  induced  by  the  process.  As 
for  the  computations  mentioned  above,  at  present  there  are  no  adequate  experi¬ 
mental  data  or  accurate  solutions  for  comparison.  Figure  6  gives  the  positions 
of  the  minimum  pressure  points;  by  comparing  the  computation  from  the  figure 
and  the  numerical  solution  of  N-S  equation  in  manuscript  [11],  we  see  that  these 
two  are  very  close.  The  drag  variations  of  pressure  difference  around  the 
cylinder  are  listed  in  Fig.  7. 
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Fig.  7.  Drag  coefficient:  _  computations  in 

the  paper;  +  Rel.535*104  experimental  value  [36]. 

O  Re  2.030x10*.  [36],  e.  Re  3.200X10*  [36J 

v  Re  3.280  X  10*.  [36J.  Q  Re  2.070X10*  [36]  - 

Comparing  the  section  (A)  computations  with  the  experiment  in  manuscript  [36], 
there  is  a  relatively  good  agreement.  Curve  B  is  slightly  lower  than  the  experi¬ 
mental  values;  however,  when  the  effect  of  the  discrete  vortexes  is  taken  into 
account  before  t«1.2,  the  drag  will  increase  and  equal  to  the  experimental 
result.  In  the  figures,  the  computational  results  of  the  concentrated  vortex 
model  are  given;  the  results  depart  considerably  from  the  actual  results.  Varia¬ 
tions  of  the  separation  points  at  the  initial  motion  period  are  listed  in  Fig.  8; 
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Fig.  8.  Separation  points  of  boundary  layer: 

______  solutions  in  the  paper; 

- rill  -X-X-*-(3l] 

this  is  the  result  of  coupling  computation  under  mutual  action  of  internal  and 
external -layer  flows.  With  increase  of  time,  at  first  the  separation  points 
rapidly  move  forward,  and  then  vary  slowly.  Computations  in  the  paper  show 
results  in  the  initial  stage  consistent  with  those  in  [28];  however,  the  method 
is  relatively  simple  and  relatively  close  to  the  numerical  solution  of  N-S  equation 
in  [11].  Comparing  computations  in  the  paper  and  the  Polhausen  method  in  [31], 
there  is  an  obvious  improvement;  the  divergent  angle  computed  from  the  Polhausen 
method  is  obviously  understated.  At  t«1.2,  the  differences  between  the  separation 
points,  and  the  studies  [28]  and  [31]  are  due  to  different  separation  criteria. 
Curve  B  differs  only  slightly  from  curve  C;  this  means  that  the  effect  of  the 
discrete  vortex  produced  when  t<1.2  is  mainly  limited  to  the  neighborhood  of  the 
rear  stationary  point,  having  little  effect  on  most  areas  around  the  cylinder. 

This  is  consistent  with  those  shown  in  Fig.  5,  the  flow  characteristics  at  an 
object  surface. 

This  paper  is  limited  to  studies  of  the  initial -period  situation  of  separation 
flow;  the  motion  time  corresponding  to  the  computations  is  relatively  short.  If 
computations  are  continued  with  increase  in  time,  the  separation  points  will 
slowly  move  forward  (Fig.  8);  the  drag  also  increases  with  the  continuous  enlarge¬ 
ment  of  the  symmetric  vortexes  (as  shown  by  curves  B  and  C  in  Fig.  7).  This  trend 
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is  not  consistent  with  actuality.  Actually,  after  the  motion  time  is  increased, 
the  shearing  layer  behind  the  cylinder  will  be  separated  and  reverse-direction 
vortexes  are  generated  under  the  action  of  reverse-pressure  gradient.  Because 
the  buildup  and  influence  of  the  reverse-direction  vortex  have  a  bearing  on 
distributions  of  pressure  and  velocity,  as  well  as  on  the  separation  point  and 
drag,  the  above-mentioned  computations  in  the  paper,  thus  far,  do  not  include 
this  factor,  which  will  be  discussed  in  another  paper. 

Conclusions 

This  paper  uses  a  method  of  combining  the  discrete  vortex  model  and  boundary 
layer  theory  to  study  the  unsteady  initial-stage  separation  flow  in  a  situation 
of  high-Reynolds-number,  incompressible  laminar  flow.  The  study  derives  the 
variation  rule  Cwith  time)  of  vortex  motion,  distribution  of  pressure  and  velocity 
of  flow  field,  and  separation  points  of  boundary  layer  and  drag.  The  vortex 
buildup  process  obtained  by  computations  closely  resembles  the  diagrams  plotted 
experimentally.  Variations  of  the  distribution  of  pressure  and  velocity  of  the 
object  surface  are  rational;  compared  to  available  results,  computations  of 
minimum  pressure  points  are  very  close.  In  coupling  computations  of  separation 
points,  the  Stratford  method  is  applied  to  the  quasi-steady  boundary  layer  of 
external  flow  as  varying  with  time;  the  computations  are  simple,  yielding  satis¬ 
factory  results.  The  drag  computations  are  quite  consistent  with  the  experimental 
results. 

The  detailed  investigations  in  the  paper  on  flow  characteristics  are 
advantageous  to  the  analysis  of  flow  regime.  The  variations  of  distributions 
of  pressure  and  velocity  at  the  object  surface  reveal  that  after  a  short  time 
following  the  circulation  flow  zone,  the  shearing  layer  of  the  motion  is  under  the 
control  of  the  reverse-pressure  gradient;  attention  should  be  given  to  the 
separation  of  the  rear  shearing  layer  and  its  effect  on  flow. 

Before  the  boundary  layer  begins  its  quasi-steady  variations,  the  discrete 
vortexes  produced  by  the  initial  separation  exert  effects  on  vortex  motion  behind 
the  cylinder  and  on  the  pressure  and  velocity  distributions  in  the  local  zone  at 
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the  object  surface.  Moreover,  the  discrete  vortexes  slightly  increase  the  drag 
this  has  little  effect  on  the  separation  points.  The  initial  positions  and 
shedding  times  of  these  discrete  vortexes  should  be  determined  from  the  M.R.S. 
separation  criterion. 
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At  present,  most  of  the  methods  for  solving  the  complete  three- 
dimensional  flow  in  turbomachinery  have  their  foothold  on  Wu's  theory 
of  two  families  of  stream-surfaces.  In  this  paper  the  concept  of 
the  accuracy  of  the  complete  three-dimensional  solution  is  clarified. 
The  specification  of  boundary  conditions  is  discussed  for  the  subsonic 
case  and  a  new  suggestion  for  giving  the  boundary  condition  is  pre¬ 
sented.  It  is  shown,  that  only  if  the  boundary  conditions  of  the 
two  families  of  stream  surfaces  are  fully  compatible  in  the  sense 
of  rigorously  satisfying  the  three-dimensional  equations,  would  the 
complete  three-dimensional  flow  solution  converge  to  in  an  error 
less  than  any  specified  small  value.  In  connection  with  the  analysis 
of  this  paper  some  preliminary  comparison  and  evaluation  have  been 
made  with  regard  to  the  methods  for  solving  complete  three-dimensional 
flow  in  turbomachinery  published  since  1978 . 


Foreword 


Advances  have  been  made  at  home  and  abroad  concerning  the  three-dimensional 
flow  in  turbomachinery;  the  so-called  complete  three-dimensional  solution  was 
introduced  [1,  2],  However,  this  solution  is  still  imperfect.  For  example,  in 
calculating  the  internal  flow  field  of  a  longitudinal -flow  turbine,  the  result 
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in  the  study  [1]  presents  a  considerable  warping  stream  surface,  concluding  the 

necessity  for  a  complete  three-dimensional  solution.  However,  the  results  in  the 

report  [2]  presents  very  little  warpage  of  the  s^  stream  surface  with  very  little 

difference  from  the  conventional  two-dimensional  solution.  Likewise,  the  study  [1] 

concludes  that  the  accuracy  of  velocity  convergence  is  limited  to  a  numerical 

solution;  the  accuracy  is  different  with  different  physical  problems,  as  low  as 

6  percent  in  the  examples  cited  in  the  manuscript.  However,  the  report  [2]  does 
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not  have  this  limitation,  arriving  at  an  error  of  3-10  for  the  relative  Mach 
number.  Which  one  is  correct?  With  reference  to  our  work  [3],  these  problems 
are  analyzed  here  with  relatively  detailed  comparison  and  discussion. 

Discussions  on  the  Method  of  Solution 

The  solution  methods  of  [1],  [2]  and  [3]  are  based  on  Professor  Wu  Zhonghua's 
theory  of  two  families  (s^  and  s^)  of  stream  surfaces.  This  theory  serves  to  reduce 
the  stream- surface  concept  of  a  three-dimensional  problem  to  a  two-dimensional 
problem  of  s^  and  s^  stream  surfaces.  It  can  be  proved  that  these  two-dimensional 
problems  are  described  by  nonlinear  elliptical  type  partial  differential  equations, 
under  the  premise  of  steady,  non-viscous  and  subsonic  flow.  A  complete  three- 
dimensional  solution  aims  at  finding  a  three-dimensional  field,  which  can  satisfy 
a  series  of  boundary  problems  at  s^  and  s^  stream  surfaces.  We  introduce  the 
following  descriptions  in  the  visualization:  a  certain  blade-channel  space  is 
divided  into  a  series  of  s^  stream  surfaces  with  a  certain  thickness  of  a  stream 
slice,  as  well  as  for  s^  stream  surfaces.  To  a  two-dimensional  problem,  two  geo¬ 
metric  properties  of  shape  and  thickness  are  included  in  the  solution  of  the 
boundary  problem  at  a  stationary  stream  surface.  However,  the  three-dimensional 
problem  involves  solving  a  boundary  value  problem  of  two  floating  stream  surfaces. 
The  solution  of  each  series  of  stream  surfaces  relates  to  another  series  of  stream 
surfaces  in  how  their  geometric  properties  must  be  modified.  With  continuation 
of  alternate  superposition,  the  variables  decrease  continuously  until  complete 
motionless.  In  this  case,  the  solution  of  two  sets  of  the  boundary  value  problems 
will  coincide  within  a  three-dimensional  flow  field;  this  is  the  final  solution  of 
a  complete  three-dimensional  flow.  Obviously,  it  does  not  share  the  same  concept 
as  the  conventional  two-dimensional  boundary  problem;  it  is  entirely  different  as 
to  the  method  of  determining  the  solution  accuracy. 
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According  to  the  discussion  above,  the  following  conditions  can  be  used  to 
evaluate  the  accuracy  of  a  complete  three-dimensional  solution  [1,  3]: 

Sir  <e 

In  the  equation,  and  W  ^  are,  respectively,  the  relative  velocities 
found  by  solving  the  and  s^  stream  surfaces  at  the  same  nodal  point.  The 
condition  for  evaluating  the  convergent  accuracy  adopted  in  the  report  [2]  are: 

Sir  -  \  w.+,  -  <a 

In  the  equation,  v  is  the  numbers  of  superpositionings  and  w.  —  (w,  +  w,^)J2. 

Obviously,  convergence  of  W  and  does  not  indicate  coincidence  of  the 

two  sets  (s^  and  s^)  of  the  two-dimensional  solutions.  The  velocity  field  thus 
constructed  does  not  satisfy  the  boundary  value  problems  of  either  the  s^  or  s2 
stream  surface.  The  report  [2]  considers  that  this  is  the  simultaneous  satisfac¬ 
tion  in  the  "average  sense";  however,  the  "average"  itself  means  the  artificial 
weighting  randomness.  For  example,  we  can  take  W  =A.W  +A^W  (A  and  A  are 

^  ^  Sj  <.  S2  1  2 

the  weighting  constants  taken  randomly) .  It  can  be  expected  that  according  to 

the  solution  method  in  [2],  the  result  of  the  same  "convergence"  is  entirely  dif- 
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ferent.  Therefore,  further  discussion  is  needed  on  the  accuracy  of  3' 10  given 
in  [2].  In  [2],  other  series  of  the  "averaging"  method  (there  are  continuity 
equations  in  addition  to  equations  of  motions)  are  adopted.  This  choice  may  blur 
the  "individuality"  of  the  s^  and  s9  stream  surfaces,  so  the  amount  of  warpage 
of  the  Sj  stream  surface  finally  computed  is  very  small. 

Suggestion  for  the  Boundary  Value  Problem 

In  solving  the  boundary  value  problem,  for  the  boundary  condition  and  the 
actual  physical  problem  to  be  consistent,  in  the  traditional  method  of  calculating 
the  blade  grid,  the  blade  channel  is  extended  to  the  upper  and  lower  streams  from 
approximately  twice  the  grid  distance  to  a  point  at  infinite  distance  (i.e.,  the 
station).  At  this  point,  a  homogeneous  stream  condition  (Fig.  1)  is  given  in 
order  to  automatically  satisfy  the  equation  of  motion  at  the  boundary.  Thus,  the 
physical  problem  itself  can  lead  to  a  rational  suggestion;  however,  the  shape  of 
the  boundary  stream  and  the  distribution  of  stream  slice  thickness  values  are 

23 


difficult  to  be  given  in  advance.  For  example,  in  calculation  of  stream  surface, 
to  the  upper  stream  station  of  the  axial-flow  type  compressor  with  internal  gas- 
inlet  cone  and  to  the  lower-stream  station  of  centrifugal -type  compressor,  the 
randomness  of  extension  of  the  meridian  channel  is  very  obvious.  For  a  centrifugal 
turbine,  in  order  to  avoid  supersonic  difficulties  while  increasing  the  radius, 
the  median  channel  in  the  lower  stream  must  be  changed  to  the  expansion  type. 
Besides,  in  the  intermediate  stage  of  a  turbo-machine,  as  contrasted  to  a  wind 
tunnel  with  a  plane  blade  grid,  often  another  row  of  blades  is  located  with  a 
very  small  gap  in  the  front  and  rear  of  the  other  row  of  blades.  Similarly,  the 
boundary  streamlines  are  difficult  to  extend  toward  the  upper  and  lower  streams. 

We  calculated  two  similar  centrifugal  turbines;  the  only  distinction  is  that  the 
difference  of  x-coordinates  (of  the  meridian-type  lines  outside  the  (K+l)  station) 
is  just  0.72  mm.  If  it  is  considered  that  the  (K+l)  station  is  a  lower-stream 
extension  station,  the  calculation  indicates  that  the  change  is  very  obvious  as 
caused  by  the  turbine  outlet  flow  field  (K  station).  Figure  2  shows  the  velocity 
distribution  as  the  s^  flow  surface  divided  in  the  middle. 


Fig.  1. 


From  [5],  we  introduced  another  solution  method:  the  boundary  is  limited  at 
the  front  or  rear  edge  of  the  blade  grid,  or  extended  outward  a  very  short  sta¬ 
tion;  the  parameters  at  the  boundary  and  the  internal  flow  field  are  solved  simul¬ 
taneously.  For  an  example  of  ifrO  in  a  standard  two-dimensional  boundary  problem 
(Fig.  3),  if  the  same  equidistant  steps  are  taken  in  the  x  and  y  directions,  and 

the  boundary  values  are  given  at  the  upper  and  lower  wall  surfaces  for  the 
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subscripts  K*Km  and  K  ,  then  all  the  internal  node  points  (the  upper-stream  points 
with  subscripts  j=jQ,  and  the  lower-stream  points  with  subscripts  j«j  )  in  the 
solution  region  can  be  solved,  respectively,  by  the  following  central,  forward, 
and  rearward  differential  forms: 

•*  +  —  K  +  •J'l.K-l  +  "*  0 

—  -4>)e+\.Z  —  4-  </>/..*_!  4-  4>I..K+\  «  0 

'Uim-hK  —  +  4/m.K-t  +  —  0 

Altogether,  (K  -Ko-l)x(ja-jo+l)  linear  independent  equations  are  written. 
Hence,  the  solution  includes  the  unknown  function  ip  of  the  same  number  of  points 
of  the  internal  node  points,  and  upper-  and  lower-stream  points. 


Fig.  2. 


Fig.  3. 


For  a  classic  boundary  problem,  write  (Ka-K0+l)x(jm-j0-l)  equations  according 
to  the  central  differential  form  of  the  internal  node  points.  In  addition,  the 
solution  concerning  only  the  internal  node  points  is  sought.  For  the  points  in 
the  upper  and  lower  streams,  the  boundary  values  are  given  according  to  some 
variation  rule  of  the  unknown  function  (the  first  category)  or  its  derivative 
(the  second  category).  In  comparison,  we  can  see  that  our  solution  method  only 
gives  the  boundary  condition  corresponding  to  the  variation  rule  of  gas -flow 
parameters  required  by  the  equation  of  motion.  It  should  be  emphasized  that  the 
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above-mentioned  differential  form  corresponds  to  the  terminal -point  condition  of 
the  constant  second-order  derivative  when  numerical  calculations  are  made  according 
to  the  sample-stripe  function.  As  indicated  by  computations  and  practice,  for  a 
two-dimensional  problem  the  suggestion  of  this  boundary  condition  is  still  steady 
and  convergent.  Using  a  streamline  iterative-substitution  method  on  a  computer, 
satisfactory  accuracy  as  to  the  convergence  of  the  entire-field  streamlines 
(including  the  inlet  and  exit  stations)  can  be  obtained;  for  example,  the  absolute 
error  is  less  than  10  ^  mm  and  during  the  iterative- substitution  solution,  no 
effect  on  the  perturbation  of  the  initial  value  is  shown.  We  believe  that  this 
suggestion  method  leads  to  a  relatively  rational  physical  problem  because  the 
equation  is  also  satisfied  at  the  boundary  (although  the  central  differential  form 
cannot  be  used,  the  equation  is  not  strictly  satisfied).  Moreover,  the  randomness 
of  extending  to  the  upper  and  lower  streams  can  be  avoided;  therefore,  this  method 
is  appropriate  as  a  treatment  method  for  boundary  conditions  in  engineering  calcula- 
t ions . 

This  solution  method  is  used  as  the  complete  three-dimensional  solution  in 
the  report  [3],  In  addition,  to  satisfy  the  periodic  condition  at  the  turbine 
exit,  we  extended  one  station  from  the  turbine  exit  as  the  lower-stream  boundary. 

The  amount  of  adjustment  of  the  circumferential-direction  coordinate  9  is  expressed 
by  A0;  the  difference  of  velocities  between  the  suction  and  pressure  surfaces  on 
the  same  revolving  arc  is  expressed  by  M.  These  two  quantities  are  similar  to 
those  in  [5] .  We  may  adopt  any  method  for  solving  roots  to  automatically  find 
A9,  on  a  computer,  by  satisfying  AW=0.  The  report  [3]  makes  detailed  three-dimen¬ 
sional  flow  field  computations  for  a  rearward  bent  turbine.  Figures  4  and  5  give, 
respectively,  the  ratios  of  streamline  shapes  and  velocity  distributions  before 
and  after  the  boundary  adjustment  of  the  upper  and  lower  streams  of  the  s^  stream 
surface  that  is  the  closest  to  the  external  wall.  Moreover,  we  are  obliged  to 
point  out  that  as  the  periodic  conditions  are  increasingly  satisfied,  the  compati¬ 
bility  of  boundary  values  of  the  and  s2  stream  surfaces  is  also  improved, 
thereby  raising  the  accuracy  of  the  three-dimensional  flow-field  solution  from 
6^*10  percent  before  the  adjustment  to  ^=3.5  percent. 

We  can  see  that  in  the  computations  in  [2],  no  such  adjustment  is  conducted. 

In  the  velocity  fields  given  in  Figs.  7,  8  and  9  in  manuscript  [2],  the  greater 

26 


the  loading,  the  closer  is  the  exit.  So  the  s2  stream  surface  does  not  remain 
close  to  the  pressure  surface  of  blades  as  seen  in  actual  phenomenon;  conversely, 
the  s2  stream  surface  is  close  to  the  suction  surface  of  blades  (the  sliding 
coefficient  u  of  the  turbine  will  be  greater  than  1). 

We  have  to  point  out  that  the  above-mentioned  adjustment  method  often  is 
not  effective  for  centrifugal  turbines  with  high  loading  at  blades  and  constantly 
in  separation  flow  condition  as  proven  in  experiments,  especially  for  turbines 
with  very  few  blades  at  the  diametral-direction  exit  or  with  few  blades.  As  indi¬ 
cated  by  computations,  in  this  case  the  deviation  amount  of  the  streamline  position 
at  the  lower- stream  boundary  is  so  large  that  it  leads  to  too  great  a  variation 
in  streamline  curvature  at  the  exit  station  because  of  the  non-viscous  calculation, 
that  the  blade  surface  is  required  to  be  the  stream  surface.  Thus,  computations 
often  result  in  overflow.  This  also  explains  the  relationship  between  the  physical 
phenomenon  and  the  mathematical  description;  i.e.,  separation  flows  will  exist  in 
the  neighborhood  of  the  blade  exit.  When  a  flow  function  equation  is  adopted  in 
seeking  the  solution,  at  stations  with  more  frequent  computations  of  the  tail- 
edge  region,  the  calculation  results  can  satisfy  the  periodic  condition;  however, 
forced  unloading  will  certainly  lead  to  abrupt  rises  in  pressure  (often  accompany¬ 
ing  with  violent  oscillating  motions)  on  the  blade  suction  surface  in  the  neighbor¬ 
hood  of  the  blade  exit.  Actually,  flows  have  been  separated  some  time  ago; 
therefore,  the  calculation  result  is  still  noncredible.  Examples  in  [1]  show 
smaller  turbine  loads  (the  difference  of  the  relative  Mach  numbers  of  two  blade 
surfaces  at  the  exit  cross-section);  however,  whether  or  not  the  calculations  are 
consistent  with  the  actual  model  depends  on  further  proof. 

Accuracy  of  Complete  Three-Dimensional  Solution 

Actually,  the  accuracy  of  a  complete  three-dimensional  solution  is  the  extent 
of  coincidence  of  two  families  of  the  two-dimensional  boundary  value  problem. 
Therefore,  the  first  factor  affecting  accuracy  is  the  compatibility  of  the  boundary 
condition.  Moreover,  since  and  s2  are  always  solved  separately;  their  dif¬ 
ferential  equations  do  not  approach  the  original  equations  with  the  same  charac¬ 
teristics.  For  example,  since  the  grid  width  ratios  of  s^  and  s2  directions 

change  greatly  (in  the  case  of  diametral -flow-type  turbine)  from  inlet  to  exit, 
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it  is  impossible  to  hsve  the  same  interruption  error.  This  is  the  reason  why 
the  accuracy  is  limited  to  a  certain  extent.  The  report  [1]  considers  that  the 
reason  is  the  transformation  of  local  grid;  this  is  the  effect  of  the  second 
factor.  The  authors  consider  that  the  first  factor  may  be  the  inherent  reason. 

First,  let  us  view  a  very  interesting  numerical  experiment:  from  example  2 
of  the  report  [3],  when  we  adopt  the  solution  method  in  which  both  inlet  and  exit 
stations  are  included  in  the  solution  range,  and  about  40  iterative  substitutions 
are  conducted  at  s^  and  s^,  the  "best"  accuracy  thus  obtained  is 

s*-  I  (w\  -n\) 5%. 

The  convergent  accuracies  of  the  streamline  grid  at  0  and  n  directions  are 

ec  =0.3  mm.  By  further  iterative  substitutions,  very  slow  divergence  will 
sl’s2 

occur  (refer  to  Fig.  6  for  the  computation  results  with  a  different  axial-flow 
turbine).  However,  if  we  fix  the  result  of  the  40th  iterative  substitution  accord¬ 
ing  to  the  streamline  distribution  at  the  inlet  and  exit  stations  (corresponding 
to  the  boundary  condition  of  not  strictly  satisfying  the  motion  equation)  before 
iterative  substitutions  of  and  then  the  solution  accuracies  can  be  improved 


to  eJ&  percent  and  e 

n  S  .  j  S 

iterative  substitutions 


^0.04  mm,  and  a  very  slow  increase  occurs  with  continuous 

2 

.  If  we  further  fix  the  streamlines  of  the  station  inward 


from  the  inlet  and  exit,  then  the  solution  accuracy  can  be  further  improved  to 

eJ\30.5  percent  and  e  20.03  mm  (Fig.  7).  Moreover,  in  the  computations  mentioned 
si’  2 

above,  only  when  the  maximum  error  of  the  streamline  position  occurs  at  the  boundary 
cross  section  of  the  solution  region,  is  the  "limiting"  accuracy  of  the  solution 


obtained. 


The  computations  mentioned  above  indicate  an  important  effect  of  the  boundary 
values.  The  reason  why  the  fixed  boundary  can  improve  solution  accuracy  is  because 
of  the  forced  compatibility  of  the  boundary  conditions  that  are  originally  incom¬ 
patible  for  Sj  and  $2,  and  the  flow  field  in  the  vicinity  of  the  boundary  is  made 


a)  The  authors  express  their  gratitude  to  Zhu  Jungguo  for  his  suggestions  and 
assistances  in  the  above-mentioned  computations. 
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to  adapt  to  the  boundary  condition,  in  addition  to  advantages  due  to  solving  the 
first-category  boundary  problem.  However,  since  the  boundary  conditions  are  in 
forced  compatibility,  the  improvement  of  accuracy  is  limited  because  at  that  time 
the  boundary  conditions  do  not  strictly  satisfy  the  three-dimensional  equations. 
While  the  cross-section  of  the  stationary  boundary  moves  further  inside,  the 
further  the  boundary  conditions  are  satisfied  with  regard  to  the  three-dimensional 
equations,  the  higher  is  the  compatibility.  Hence,  the  process  of  improving 
accuracy  mentioned  above  reappears  once  more  at  a  higher  level.  Therefore,  we 
consider  that  only  when  the  given  boundary  conditions  strictly  satisfy  the  three- 
dimensional  motion  equations,  are  the  conditions  completely  compatible.  At  that 
time,  the  accuracy  of  the  three-dimensional  flow-field  solution  can  be  randomly 
given.  We  also  can  deduce  the  following:  if  the  solution  region  can  be  extended 
to  an  infinite  distance  toward  the  upper  and  lower  streams,  the  flow  fields  there 
are  uniform  (in  other  words,  automatically  satisfying  the  three-dimensional  motion 
equation) ,  therefore  the  accuracy  of  the  three-dimensional  flow-field  solution  can 
be  randomly  given. 


Fig.  4.  _ after  adjustment;  Fig.  5.  _ after  adjustment; 

—  before  adjustment.  —  before  adjustment. 

We  must  explain  that  the  discussion  above  is  about  the  method  of  obtaining 
a  three-dimensional  solution  by  using  the  alternate  iterative  substitutions  of 
two  two-dimensional  problems.  This  method  is  brief  with  well-defined  concepts; 
however,  there  are  new  problems  in  the  numerical  solution.  Usually  in  a  two- 
dimensional  boundary  problem,  the  boundary  value  can  be  randomly  given  (even 
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Fig.  6  Fig.  7. 

Key:  (*)  Times  of  iterative  Key:  (*)  Times  of  iterative 

substitution.  substitution. 

though  this  does  not  necessarily  satisfy  the  equation) ,  to  arrive  at  steady  conver¬ 
gent  solutions  mathematically.  However,  in  the  three-dimensional  solution,  the 
requirement  of  compatibility  of  two  two-dimensional  problems  is  suggested;  the 
distinction  between  two-  and  three-dimensional  solutions  is  also  manifested  in  the 
effect  of  grid  density.  For  the  same  example,  along  the  streamline  direction  the 
number  of  calculation  stations  is  constant;  however,  while  solving  for  motion 
equations,  if  different  densities  of  grids  are  taken  in  the  Sj  and  s,,  directions, 
the  accuracy  is  higher  if  we  use  a  denser  grid  in  a  two-dimensional  problem.  The 
result  is  just  the  reverse  in  a  three-dimensional  problem.  Moreover,  we  know 
that  different  equations  can  have  determinate  solutions  only  when  the  initial 
and  boundary  conditions  are  given  simultaneously.  This  basic  fact  affects  the 
accuracy  of  solutions  by  mutually  iterative  substitutions  of  differential  equa¬ 
tions;  basically  this  is  different  from  the  solution  of  iterative  substitution  of 
algebraic  equations.  Previously,  we  discussed  the  accuracy  of  the  iterative 
substitution  solution  of  differential  equations  by  using  a  procedure  of  reducing 
the  number  of  dimensions.  In  recent  years,  some  researchers  used  a  procedure 
for  varying  the  number  of  unknown  functions  to  construct  iterative  substitution 
solutions  of  differential  equations.  At  that  time,  similar  limitations  also 
exist  in  the  solution  accuracy;  for  example,  in  the  solution  method  mentioned  in 
[6],  the  errors  can  be  convergent  to  less  than  6  percent. 
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Conclusions 


1.  By  using  the  alternate  iterative  substitution  for  two-dimensional  problems 
of  two  families  of  stream  surfaces,  solutions  of  three-dimensional  flow  fields  can 
be  derived.  This  solution  method  introduces  new  mathematical  problems  that  should 
be  discussed.  Our  current  understanding  is  as  follows:  the  accuracy  of  solutions 
for  complete  three-dimensional  flows  is  indicated  by  the  extent  of  coincidence  of 
solving  two-dimensional  problems  of  these  two  families  of  stream  surfaces.  This 
concept  differs  from  the  concept  of  accuracy  in  solving  the  conventional  two- 
dimensional  boundary  problems.  Only  when  the  boundary  conditions  are  completely 
compatible  while  strictly  satisfying  the  three-dimensional  equations,  can  the 
accuracies  of  solutions  of  three-dimensional  flow  fields  be  given  randomly. 

2.  By  simultaneously  solving  for  boundary  parameters  and  internal  flow  fields, 
in  other  words  by  using  the  solution  method  of  given  boundary  values  for  the 
variation  rule  required  by  the  motion  equations,  the  physical  problem  becomes 
rational,  relatively  speaking. 

3.  In  computations  neglecting  viscosity  and  separation,  starting  from  satisfy¬ 
ing  the  periodic  conditions,  it  is  necessary  to  make  adjustments  to  the  upper- 

and  lower-stream  boundaries;  this  is  also  advantageous  in  improving  the  solution 
accuracy  of  three-dimensional  problems.  However,  for  centrifugal  turbines  with 
high  blade  loads,  a  calculation  method  should  be  needed  in  a  more  realistic  model 
of  separate  flow. 
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Paper  presents  a  numerical  method  for  computing  the  various 
compressible  steady  laminar  flows  in  turbomachinery.  A  numerical 
solution  is  obtained  by  iterative  computation  based  on  the 
enthalpy  gradient  equation,  energy  equation  and  entropy  equation, 
expressed  along  arbitrary  curved  lines.  This  paper  at  first 
analyzes  and  discusses  the  computational  method  in  a  stationary 
coordinate  system,  and  describes  the  method  of  axisymmetrical 
flow  in  the  stator.  As  revealed  by  the  computational  examples, 
the  characteristics  of  the  method  are  brief  procedures  and  high 
computation  speed;  therefore,  the  method  may  be  widely  applied 
in  engineering  design. 


I .  Foreword 

Since  Professor  Wu  Zhonghua  presented  the  theory  of  stream  surfaces  [1] 
in  1952,  the  computational  method  of  three-dimensional  flow  in  turbomachinery  has 
made  rapid  strides.  At  the  present,  efforts  are  underway  to  obtain  the  real  three- 
dimensional  flow  solutions  (for  example,  in  [2]).  However,  all  these  computational 
methods  consider  fluid  in  flow  channels  as  a  nonviscous  ideal  fluid,  which  encoun¬ 
ters  mechanical  loss  due  to  its  viscosity  in  the  flow  process;  this  loss  can  be 
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compensated  by  a  loss  factor  obtained  experimentally.  Although  this  procedure 
avoids  difficulties  in  calculating  viscous  flow,  yet  certain  errors  are  introduced. 
In  addition,  this  method  cannot  be  used  to  compute  problems  of  secondary  flows  in 
turbomachinery . 

Later,  Professor  Wu  derived  a  fundamental  set  of  equations  (in  manuscript  [3]) 
that  a  viscous  compressible  flow  in  turbomachinery  must  satisfy.  In  addition,  he 
conducted  detailed  analyses  of  the  function  of  viscous  force  and  the  physical 
significance  of  the  viscous  items.  For  other  researchers,  Patankar  and  Spalding 
[4]  proposed  adopting  a  propulsion  integration  method  to  calculate  viscous  flow 
in  piping  channels;  Pratap  and  Majundar  [S]  adopted  this  method  to  calculate  the 
viscous  incompressible  flow  in  a  rectangle-shaped  cross-section  straight -pipe  chan¬ 
nel  as  the  flow  revolves  around  the  axis  of  main-flow  direction;  and  Ghia  et  al 
(6]  adopted  a  polar  coordinate  system  to  calculate  the  viscous  incompressible  flow 
in  a  sector-shaped  cross-section  of  a  straight-pipe  channel.  These  two  above- 
mentioned  pipe  channels  closely  resemble  turbomachinery  channels;  however,  limited 
by  the  forms  of  equations  they  adopted  and  because  during  calculations,  they 
neglected  the  second-order  derivative  items  along  the  main-flow  direction,  the 
method  only  solves  the  problem  of  viscous  flow  in  regular-shaped  cross-section 
pipe  channels.  Until  now,  there  has  still  been  no  solution  in  the  calculation  of 
viscous  compressible  flow  in  irregular- shaped  variable  cross-section  bent  pipe 
channels  like  turbomachinery  channels. 

This  paper  proposes  a  numerical  computation  method  for  viscous  compressible 
steady  laminar  flow  of  an  internal  flow  field.  First,  starting  from  the  motion 
equation  (Navier-Stokes  equation)  of  viscous  fluid,  am  enthalpy  gradient  equation 
is  derived.  Together  with  the  energy  equation  and  the  entropy  equation,  the 
enthalpy  gradient  equation  is  used  in  making  iterative  computations.  Through  flux 
calibration,  numerical  solutions  of  viscous  flow  of  the  entire  flow  field  can  be 
obtained.  As  equations  adopted  in  the  paper  are  written  in  the  derivative  form 
along  directions  of  three  random  curves,  the  method  can  be  used  to  compute  any 
shape  bent  pipe  channels  (including  turbomachinery  channels)  of  variable  cross 
sections  in  solving  the  viscous  compressible  steady  flow  problems. 
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At  first,  the  paper  analyzes  and  discusses,  with  the  assumption  of  a 
stationary  system,  to  derive  forms  of  equations;  these  equations  are  used  to 
compute  the  viscous  compressible  symmetric  flow  of  stationary  blade  rows  in  turbo- 
machinery.  The  situation  in  these  computations  corresponds  to  considering  the 
stationary  blade  rows  as  the  number  of  blades  approaches  infinity  and  the  blades 
become  infinitesimally  thin  to  be  an  extreme  case.  In  this  case,  the  three- 
dimensional  flow  field  of  the  stationary  blade  rows  can  be  considered  as  composed 
of  an  infinite  number  of  S2  flow  surfaces,  such  as  arc  surfaces  of  blades.  There¬ 
fore,  there  are  the  same  distributions  of  flow  parameters  on  each  stream  surface. 
In  mathematical  expression,  this  is  an  axisymmetric  flow  problem  because  flow  is 
not  related  to  0;  in  mechanics,  this  can  be  approximately  considered  as  making 
computations  on  the  intermediate  S,,  stream  surface.  In  design  problems  (reverse 
problems),  this  flow  model  is  continually  adopted  in  making  design  calculations 
of  blades.  In  checking  problems  (positive  problems),  on  this  basis  the  method  can 
be  used  to  conduct  computations  on  the  stream  surface  in  order  to  obtain  the 
numerical  solutions  of  quasi-three-dimensional  flow  in  stationary  blade  rows. 

II.  Fundamental  Equations  in  Viscous  Compressible  Flow 

1.  The  fluid  is  a  compressible  Newtonian  fluid  and  satisfies  Fourier's  heat 
conduction  flow  law.  2.  The  fluid  is  an  ideal  gas.  3.  The  flow  condition  is 
steady.  Therefore,  the  fundamental  equations  that  a  viscous  compressible  steady 
flow  should  satisfy  are 


V  •  (pv)  -  0 

O) 

•  F  —  1  v?  -+■  -  r  •  (r ! 

(2> 

P  P 

^  •  g  •+■  —  o  •  Vp  -t-  i  |tJ  •  Vu 

(3> 

P  P 

p  —  RpT 

(4> 

h  -  c,r 

(5/ 

In  the  equations,  the  unknowns  are  flow  velocity  v  of  the  fluid,  fluid  density  p, 
fluid  pressure  p,  fluid  temperature  T,  and  enthalpy  h  of  a  unit  mass  of  fluid. 
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The  mass  force  F  and  amount  of  heat  absorption  Q  are  given;  generally,  these  two 
quantities  can  be  neglected.  Also  given  are  the  gas  constant  R  and  the  constant- 
pressure  specific  heat  Cp.  Under  the  hypothesis  1,  the  viscous  stress  tensor  [t] 
and  amount  of  heat  conduction  q  satisfy,  respectively. 

Lx j  -  n  jvu  +  (Vv)7  —  1  V  •  v[I]  |  (t>) 

T  Q - k?T  (7) 

In  the  equations,  (Vv)  is  the  transposed  tensor  of  tensor  (Vv),  [I]  is  the  unit 

tensor,  y  is  the  dynamic  viscous  coefficient  of  the  fluid,  and  k  is  the  heat 
conduction  coefficient  of  the  fluid.  When  the  temperature  variation  is  small,  u 
and  k  can  be  approximately  considered  as  constants.  However,  Cp,  U  and  k  also  can 
be  considered  as  known  functions  of  temperature.  By  so  doing,  no  difficulties  will 
be  incurred  in  the  computation  though  the  computation  work  will  be  increased.  From 
the  foregoing,  equations  (1)  through  (5]  are  closed;  therefore,  theoretically 
speaking,  solutions  can  be  obtained.  However,  in  general  situations  it  is  quite 
difficult  to  solve  these  equations.  Hence,  before  we  conduct  numerical  computa¬ 
tions  these  equations  are  first  transformed. 

In  equation  (1),  its  integration  form  can  be  used  in  this  method. 

|  j  pv  •  dA  «•  JO 

In  the  equation,  A  is  the  area  of  any  cross  section  in  the  flow  channel  and  Gq  is 
the  given  flow  quantity. 


To  equation  (2) ,  by  using  the  fundamental  thermodynamic  equation 


TdS  —  dh-  -  dp 
P 


or  its  gradient  form 

TVS  —  vh  —  -Vp.  (10) 

P 

Combining  equations  (2)  and  (10),  cancelling  out  the  item  ~VP  and  neglecting 

the  mass  force  F,  we  derive 

Vh  -  TVS  +  -  v  •  [r]  —  (u  ■  V)o  (11) 

P 

This  equation  indicates  the  relationship  equation  that  enthalpy  gradient  should 

satisfy  in  a  flow  field;  therefore,  it  is  called  enthalpy  gradient  equation,  whose 

derivation  is  a  feature  of  this  method.  One  of  the  advantages  is  that  the  enthalpy 

gradient  equation  together  with  energy  and  entropy  equations  can  be  used  to  conduct 

alternate  super imposition  computation;  S  is  entropy  of  the  fluid  and  can  be  derived 

* 

from  integration  of  equation  (9)  ^  ^ lo  T *~*  (12) 


In  the  equation,  x  is  the  adiabatic  index  of  the  fluid. 

To  equation  (3),  the  equation  (10)  is  projected  along  the  streamline  direction, 
we  get 


y  •  TVS  —  v  ■  VA - -  v  •  Vp 

P 

Substituting  equation  (3)  into  the  above  equation,  we  get 

(13) 

(O') 


■  TVS  —  Q  4-  i  V‘T  +  -  [r]  • 


Vv 
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In  the  equation, %  represents  the  distance  coordinate  along  the  spatial  streamline, 
dS/d£  represents  the  variation  rate  of  entropy  along  the  streamline  £.  Therefore, 
the  equation  is  called  the  entropy  equation. 

Finally,  look  at  the  form  of  the  enthalpy  gradient  equation  (11)  along  the 
streamling  direction 

o  •  VA  —  u  •  TVS  +  —  w  •  (V  •  (.*])  —  »•(»•  V)i> 

P 

By  using  equation  (13)  and  vector  computations,  we  get 

v  •  V  (a  +  -N)  -  Q  +  -  V‘T  +  -V  •  ([r]  •  u).  (H) 

\  Z/P  P 

2 

By  introducing  ho=[h+(v  /s)]  which  represents  the  total  enthalpy  (total  energy) 
of  a  unit  mass  of  fluid,  equation  (14)  can  be  rewritten  as 

4a-±\t  +  l,.T  +  ±,.u,}-c)\  (»•> 

This  equation  indicates  the  variation  rate  of  total  energy  of  the  fluid  along  a 
spatial  streamline  Z,  therefore  equation  (14*)  is  called  the  energy  equation.  We 
can  see  from  the  above  that  the  energy  equation  is  only  the  expression  form  of 
enthalpy  gradient  equation  along  the  streamline  direction. 


37 


Therefore,  after  transformation  new  equations  are  composed  of  equations  (4) , 
(S),  (8),  (11),  (12),  (13),  and  (14);  equation  (11)  is  a  form  resulting  when  only 
two  directions  are  taken,  except  for  the  streamline  directions.  In  the  equations, 
the  other  unknowns  are  v,  p,  P,  T,  h,  and  S;  therefore,  the  equations  can  be  still 
considered  as  closed;  they  can  be  used  to  conduct  numerical  computation. 

III.  Equation  Forms  and  Boundary  Conditions  in  Turbomachinery 

According  to  the  flow  characteristics  of  fluid  in  turbomachinery,  a  discussion 
can  proceed  by  selecting  a  cylindrical  surface  coordinate  system  (r,6,  2);  the  z 
axis  is  coincident  with  the  blade  axis.  To  the  axisymmetric  flow  considered  in 
the  paper,  the  partial  derivative  of  all  flow  parameters  to  0  is  0.  Therefore, 
the  above-mentioned  equations  can  be  rewritten  as  the  following  forms: 

(A)  The  entropy  equation  (11)  can  be  rewritten  as 


dr  dr  P  I  3  dr1  3  r  dr  3  r2 

dlu,  1  I  dv,  Vg  j 

dz1  +  3  drdz  j  ~  Vm  dm  ~  r  j 
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dz  T  dz  +  p  [  dr1  "**  r  dr  3  dz3  3  t  dz 

+  i  J!sd  _  r- 
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(H)* 


In  the  equations,  d/dm  stands  for  the  total  derivative  of  flow  parameters  to  the 
meridian  streamline  when  the  parameters  vary  along  spatial  streamlines;  vr,  vQ , 
and  vz  are  components  of  fluid  velocity  v  in  the  cylindrical  coordinate  system 


*  For  axisymmetric  flow,  Therefore,  the  equations  can  satisfy  the  axisym 

metric  flow  in  a  blade-less  flow  channel.  As  for  the  blade-row  channel  discussed 
at  present,  the  item  -Jj-  should  be  retained  for  the  time  being  because  of  the 
existence  of  blade  force. 
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(r,  0,  z);  and  is  the  projection  of  v  on  meridian  plane.  The  following  rela¬ 
tionships  are  satisfied  among  v  ,  v„  and  v  : 

r  o  z 

v,  —  c^sina  v,  —  t^cosa  u„  u  sin  $  vm  — *  u  cot/3 

In  the  equations,  a  is  the  meridian  flow  angle  and  8  is  the  circumferential  flow 
angle  (see  Fig.  1) . 

For  any  spatial  curve  q  in  the  flow  field,  the  direction  derivative  of 
enthalpy  h  along  curve  q  is 
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In  the  equation,  q  can  be  arbitrarily  assumed  as  three  random  directions  without 
a  common  plane  in  order  to  derive  the  equation  satisfying  the  three  directions. 
Usually,  these  three  directions  can  be  selected:  first,  it  is  the  streamline 
direction;  the  enthalpy  gradient  equation  along  this  direction  is  the  energy 
equation  (14)  derived  previously.  Second,  it  is  the  9  direction.  For  axisym- 
metric  flow,  the  equation  along  the  9  direction  is  not  necessary.  Third,  it  is 
the  spatial  curve  without  a  common  plane  of  the  streamline  and  the  9  direction. 
This  curve  is  called  the  computation-station  line,  briefly  called  the  station 
line.  When  we  compute  on  a  middle  S£  stream  surface,  this  station  line  must 
be  situated  on  a  middle  stream  surface;  the  projection  curve  of  the  station 
line  on  the  meridian  plane  is  called  a  meridian  station  line.  Generally,  a 
meridian  station  line  can  be  taken  as  a  straight  line  approximately  perpendicular 
to  the  meridian  streamline;  this  is  the  so-called  quasi-normal  line  (Figs.  2,  3) 
in  [7],  In  this  case,  equation  (16)  can  be  rewritten  as 
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In  the  equation,  d/dn  stands  for  the  total  derivative  of  flow  parameters  to  the 
meridian  station  line  n  when  the  flow  parameters  vary  along  the  spatial  station 
line;  <f>  is  the  angle  of  the  meridian  station  line. 


(B)  The  energy  equation  (14)  can  be  rewritten  as 
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(18) 

(C)  The  entropy  equation  can  be  rewritten  as 
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All  the  three  above-mentioned  equations  are  rewritten  as  the  derivative  form 
of  m  and  n  along  directions  of  random  curves.  If  the  right-hand  items  of  the 
equations  are  predetermined,  then  these  equations  become  simple  differential 
equations.  The  problems  become  very  simple;  it  is  easy  to  make  iterative  computa¬ 
tion.  This  is  a  characteristic  of  the  method. 

In  order  to  compute  the  viscous  compressible  axisymmetric  flow  problems  at 
stationary  blade  rows  in  turbomachinery,  the  following  boundary  conditions  are 
required: 

(1)  In  the  blade-row  flow  channel,  the  velocity  at  the  internal  and  external 
wall  surfaces  must  be  0;  i.e. ,  vwall=0-  0°  the  wall  surface,  there  are  two 
temperature  situations:  1.  the  isothermal  wall  surface,  i.e.,  given  the  temperature 
distribution  situation  on  the  wall  surface;  and  2.  the  adiabatic  wall  surface,  i.e., 
the  normal -direct ion  derivative  on  the  wall  surface  should  be  0.  According  to  the 
relationship  equation  (5)  that  satisfies  the  relationship  between  temperature  and 
enthalpy,  the  boundary  conditions  of  enthalpy  can  be  suitably  written  out. 

(2)  When  calculating  flow  in  turbomachinery  blade  rows,  in  order  to  eliminate 
the  effect  and  errors  introduced  by  randomly  selecting  the  boundary  conditions 

at  the  exit,  it  is  necessary  to  extend  a  distance  toward  the  upper  and  lower 
streams  from  the  inlet  and  outlet.  In  the  region  extending  from  the  upper  and 
lower  streams,  there  is  no  function  of  blade  force  since  this  is  the  blade-less 
region.  Therefore,  in  the  middle  equation  of  equations  (15),  let  Qh-  —  o 
in  order  to  get 

djllA  _  4-  1  -  v*  4.  W±.)  (20) 

dm  pVm  V  dr‘  r  dr  r 1  dz’  > 

This  equation  satisfies  the  circumferential -direction  velocity  V  of  gas  flow  in 
the  extending  region. 
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Fig.  1. 


Fig.  2. 

Key:  (1)  Meridian  streamline; 
(2)  Meridian  station  line. 


Key:  (1)  Meridian  station 

line;  (2)  Stationary  blade 
row;  (3)  Meridian  streamline. 


IV.  Numerical  Computation  Method  and  Computation  Examples 

This  paper  adopts  the  streamline  iterative  method  for  the  numerical  computa¬ 
tion.  By  taking  examples  of  axisymmetric  flow  in  the  stationary  blade  rows  in 
turbomachinery,  the  computation  steps  in  the  method  are  described  as  follows: 

(1)  At  the  middle  S2  stream  surface  (including  the  internal  portion  of 
stationary  blade  rows  and  the  extending  region  toward  the  upper  and  lower  streams), 
first  several  station  lines  q  are  selected.  Then,  according  to  some  distribution 
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principle,  assumed  streamlines  l  are  marked.  Thus,  streamlines  and  station 
lines  constitute  a  computation  grid;  on  the  grid  points,  assume  the  distribution 
situation  of  V  ,  p ,  and  T. 


(2)  The  computation  grid  at  the  middle  S2  stream  surface  is  projected  onto 
a  meridian  plane  to  get  a  meridian  gTid  composed  of  meridian  streamline  m  and 
meridian  station  line  n.  Using  numerical  differentiation,  we  get  the  meridian 
flow  angle  a  and  meridian  station  line  angle  $ : 

a  ■"  arc  tg  I 

d*  >  along  m 


d* 

•ire  tg -  , 

dr  along  n 


(3)  In  the  internal  region  of  the  stationary  blade  rows,  according  to  the 
shape  of  arc  surface  of  the  stationary  blade  rows,  we  then  calculate  Vg  and  the 
circumferential -flow  angle  8 : 


V„  -  IV— I 

dm  1  along  streamline  1 

9  —  “"g  ~~ 

*  m 


(21) 

(22) 


For  the  reverse  problem,  Vg  can  be  computed  from  the  selected  stream  types.  Then 
the  value  of  8  can  be  computed  from  equation  (22) . 


In  the  extending  region,  numerical  integration  and  iterative  computation  on 

equation  (20)  can  be  conducted  to  obtain  the  value  of  V  .  From  equation  (22),  the 

0 

value  of  8  can  be  obtained. 


(4)  By  using  numerical  differentiation,  obtain  the  values  of  derivatives  of 
flow  parameters  V  ,  V  ,  V  and  T  along  the  directions  of  meridian  streamline  m 
and  meridian  station  line  n.  Then,  use  the  characteristics  of  derivatives 


d(_  _ 

a/ 

+ 

d* 

dm  • 

dr 

dm 

dz 

dm 

dl 

3/ 

dr 

di 

d* 

dn 

dr 

In  + 

dz 

dn 

(23) 
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In  the  equations,  £  stands  for  a  certain  flow  parameter. 


dr  — 

sin  a 

da  ^ 

COSO 

dm 

dm 

dr  ^ 

cos  q? 

da  — 

sin  <p 

dn 

dn 

These  are  known  quantities.  Then  the  values  of  3f/3r  and  3f/3z  can  be  obtained 
by  solving  the  linear  algebraic  equations  (23). 

For  second-order  derivative  values,  similarly,  by  solving  the  linear 
algebraic  equations 


£l  _ 

y+ 

(J3..+ 

d>f  ' 

)£. 

JfT 

+  %-(. 

dm1 

dr 1  V  dm) 

V  drdz 

dzdr  > 

'  dm 

dz1  V 

df  d‘r 

df  d*z 

*  dr  dm 1  * 

dz  dm1 

•• 

2U4-Y  + 

[-2L  + 

J1L) 

| 

+ 

dn * 

\  dn  / 

'  drdz 

dzdr  ) 

dn 

df  d‘r 

df  d‘z 

f  9r  d*1  + 

dz  dn 1 

JL L 

—  d*f  dr  dr 

+  JH- 

dr  dz 

•4-  - 

da  dr 

dmdn 

dr1  dm  dn 

drdz 

dm  dn 

dzdr 

dm  dn 

d1)  dz 

da  df 

dlr 

d*s 

dz1  dm  dn  dr  dmdn  dz  dmdn 
<Pf  _  »/  dr  dr  +  dr_  dx_  +  9lf  dz  dr_ 

dndm  dr1  dn  dm  drdz  dn  dm  dzdr  dn  dm 

d*f  dz  dz  df  <Pr  df  d‘z 

dz1  dn  dm  dr  dndm  dz  dndm 


(24) 


m  m  dh 

to  obtain  the  second-order  derivative  values  — , 

of  the  flow  parameters  because  theoretically,  there  must  be 

d‘i  _  a>f 

drdz  dzdr 

Therefore,  we  can  take  the  average  value  as  the  value  of 
various  equations.  44 


and 


d>f 

dzdr 


-&L 

drdz 


to  be  used  in 


i 


(5)  At  every  streamline  (except  the  wall-surface  streamlines),  perform 
numerical  integration  on  energy  equation  (18)  and  entropy  equation  (19).  Moreover, 
determine  the  integration  constant  by  using  the  values  of  total  enthalpy  and 
entropy  at  the  inlet.  Thus,  the  distribution  of  total  enthalpy  and  entropy  can 

be  obtained  for  the  entire  flow  field  (except  the  wall  surface). 

(6)  At  each  station  line,  perform  numerical  integration  on  the  enthalpy  gra¬ 
dient  equation  (17).  In  addition,  using  the  boundary  conditions  of  wall  surface, 
the  distribution  of  enthalpy  h  can  be  obtained.  Then,  using  definition  of  total 
enthalpy  and  equations  (12)  and  (4),  we  can  get  the  values  of  v,  p  and  p.  An 
integration  constant  in  numerical  integration  can  be  determined  by  using  flow 
equation  (8) . 


Fig.  4. 


(7)  According  to  the  flow  distribution  principle  determined  in  advance,  using 
new  streamline  positions  through  reverse  interpolation  of  flow  and  the  relaxation 
factor  of  less  than  1,  we  can  get  the  new  assumed  positions  of  streamlines.  Next, 
return  to  step  2  to  make  an  iterative  computation  of  the  next  cycle  until  the 
distance  between  the  assumed  position  of  the  streamline  and  new  streamline  position 
can  satisfy  the  requirement  on  accuracy.  In  this  case,  the  computation  result 
obtained  is  the  numerical  solution  that  we  require. 

In  the  above-mentioned  computation  steps,  no  detailed  problems  are  mentioned 
regarding  some  numerical  computations  and  the  writing  of  programs.  These  details 
will  be  described  and  discussed  thoroughly  in  another  paper. 
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By  using  this  method,  computations  of  several  flow  channel  arrangements 
were  conducted.  Here  only  a  part  of  the  computational  results  of  stationary 
blade  rows  of  a  certain  axial -flow-type  turbine  are  plotted  in  Fig.  5.  From  the 
figure,  we  can  see  the  effects  on  velocity  and  temperature  resulting  from 
viscosity. 


Fig.  5(a).  Distributions  of  meridian 
streamlines  and  meridian  station  lines 
in  stationary  blade-row  channel  of 
axial -flow-type  turbine. 

Key:  (a)  Inlet  extension  region;  (b) 
Stationary  blade-row  region;  (c)  Exit 
extension  region. 


Fig.  5(b).  Distributions  of  Fig.  5(c).  Distributions  of 

M  and  T  at  inlet  cross-section.  M  and  T  at  exit  cross-section. 


The  characteristics  of  this  method  are  simplicity,  intelligibility,  brief 
computation,  and  fast  computation  speed.  From  the  above  computation  examples, 

21  streamlines  and  14  station  lines  are  taken.  Assuming  the  requirement  that  the 
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streamline  error  does  not  exceed  0.1  mm,  the  computation  time  on  a  719  computer 
(125,000  operations  a  second)  does  not  exceed  five  minutes  (excluding  the 
printout  time).  Moreover,  the  internal  memory  of  the  computer  required  does 
not  have  to  be  large.  Therefore,  this  method  can  be  widely  applied  in  engineer¬ 
ing  design. 

V.  Conclusions  ; 

The  paper  proposes  a  new  numerical  computation  method.  Through  iterative 
computations  among  enthalpy  gradient  equation,  energy  equation  and  entropy 
equation  in  the  direction  of  a  random  curve,  we  can  get  numerical  solutions  of 
the  viscous  compressible  steady  laminar-flow  problem  of  any  random  shape  of 
variable  cross-section  bent  pipe  channel.  Characteristics  of  the  method  include 
the  establishment  of  fundamental  equations  of  viscous  flow  in  a  grid  of  random 
curves.  In  addition,  in  iterative  computations,  these  equations  can  be  simplified 
into  conventional  differential  equations,  thus  simplifying  the  computations. 

The  paper  discusses  the  fundamental  equations  of  a  stationary  coordinate 
system  and  the  computation  methods  of  axisymmetric  flow  problem.  Moreover,  numer¬ 
ical  computations  are-conducted  on  axisymmetric  flow  in  stationary  blade  rows 
in  turbomachinery.  The  computational  result:  indicate  that  the  method  is  simple 
and  feasible.  In  Part  II,  fundamental  equations  of  viscous  compressible  flow  of 
rotational  coordinate  system  will  be  derived.  In  addition,  numerical  computations 
of  axisymmetric  flow  of  rotational  blade  rows  in  turbomachinery  will  be  conducted. 
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A  system  of  high-speed  interferograph  for  investigating 
the  interaction  of  TEA  C02  pulsed  laser  with  matter  and  other 
fast  flow  phenomena  has  been  constructed  and  applied  success¬ 
fully.  Problems  of  the  time  synchronization  and  the  optical 
match  between  high-speed  camera  and  interferometer  have  been 
solved.  We  operated  inversely  an  ordinary  CW  He-Ne  laser  to 
obtain  a  time-controllable  pulsed  light  source.  Several 
series  of  high-speed  (microseconds)  interference  photographs 
of  laser-matter  interaction  and  the  flow-field  in  supersonic 
wind-tunnel  have  been  obtained.  This  technique  may  be  used 
in  other  fast  phenomena  and  flow-field  visualization. 


FOREWORD 


In  studying  interactions  of  high-energy  lasers  and  matter,  complex  physical 
processes  occur.  Appropriate  means  are  required  to  observe  and  analyze  these 
processes  so  many  investigation  techniques  have  been  developed,  from  the  areas  of 
optics,  electricity,  spectrum,  and  X-rays  [1,  2,  3].  The  optical  technique  is  the 
most  widely  used,  because  the  technique  involves  one  of  the  most  direct  observa¬ 
tions  of  all  approaches:  it  has  been  much  refined  in  the  stydy  of  phenomena  in 
picoseconds,  nanoseconds,  and  milliseconds. 
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To  study  the  interaction  of  a  10. 6-micron-wavelength  pulsed  laser  and  matter, 
we  need  to  develop  an  optical  method  to  investigate  phenomena  during  this  inter¬ 
action.  From  the  published  reports,  one-  and  two-dimensional  studies  in  most 
experiments  are  conducted  by  simultaneous  illumination  and  plasma  growth  with  a 
stripe  or  frame-by-frame  high-speed  camera.  In  recent  years,  image-changing  tube 
high-speed  cameras  (such  as  the  Imacon  series  [4])  were  used  in  picosecond  and 
nanosecond  studies  on  light  sources  of  different  wavelengths  and  relatively  weak 
luminosity  phenomena;  the  studies  included  stripe  and  frame-by-frame  photographic 
recording.  However,  better  recordings  have  not  been  obtained  for  phenomena  of 
non-simultaneous  illumination,  weak  luminosity,  and  disturbance  of  the  medium  near 
the  interaction,  such  as  air  shock-wave  due  to  the  interaction,  and  intensive 
vaporization  following  laser  irradiation  on  a  target.  In  particular,  few  recordings 
of  high-speed  interference  photographs  are  conducted  with  frame-by-frame  continuity. 
Since  high-speed  interference  photographs  with  frame-by-frame  continuity  makes 
possible  the  qualitative  and  quantitative  analysis  in  time  and  space,  we  developed 
a  high-speed  interferograph  device. 

During  its  development,  a  continuous -wave  1  milliwatt  He-Ne  laser  was  used  as 
the  light  source  of  an  interference  instrument,  and  the  interaction  of  a  dual- 
discharge  CO,  pulsed  laser  with  matter  was  studied.  A  controllable  type  revolving 
lens  scanning  frame-by-frame  superhigh- speed  camera  was  used  in  making  the  photo¬ 
graphs;  a  pulsed  light  source  synchronized  with  the  camera  was  also  required. 
Although  there  are  advantages  with  the  He-Ne  laser,  such  as  good  monochromaticity, 
small  and  flexible  adjustment,  and  appropriate  luminosity,  there  is  still  no  pulsed 
type  operation  available  for  the  He-Ne  laser  [1].  We  applied  its  advantages  and 
operational  properties  in  suggesting  a  reverse-program  operating  plan  to  success¬ 
fully  modify  the  laser  into  a  controllable  pulsed  type  light  source.  We  perfected 
this  technique  and  obtained  many  highspeed  interference  photographs. 

Technical  Details  of  the  Development 

With  a  TEA  C02  pulsed  laser  [5]  with  an  output  wavelength  of  10.6  microns  and 

approximately  4  joules  of  energy  at  atmospheric  pressure,  the  output  pulse  width 

can  be  adjusted  within  a  range  of  1-100  microseconds.  Different  power  densities 

of  focusing  light  beams  can  be  obtained  with  different  lens  focal  lengths.  The 
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laser  can  be  used  in  studying  physical  phenomena  occurring  in  laser-irradiated 
matter.  Then  a  high-speed  interference  investigation  system  is  used  to  record 
the  photographs.  See  Fig.  1  for  experimental  layout 


Fig.  1.  Experimental  layout. 

Key:  (1)  Plane  crystal;  (2)  Target;  (3) 

Superhigh  speed  camera;  (4)  Matching 
objective  lens;  (5)  Light  beam;  (6)  Paral¬ 
lel  light;  (7)  He-Ne  laser;  (8)  C02  pulsed 
laser;  (9)  Power  source  of  CO.,  pulsed  laser; 
(10)  Single  steady  delay  circuit;  (11)  Pulsed 
power  source  of  Ne-He  laser;  (12)  Control 
platform  of  superhigh  speed  camera. 


1.  Fast  processes  are  studied  in  the  microsecond  magnitude.  For  an  accurate 
determination,  first  the  phenomena  in  the  study  must  be  fully  recorded  within  a 
short  duration.  All  investigation  segments  should  be  time-synchronized.  The 
shorter  the  duration,  the  more  the  synchronized  segments,  thus  introducing  more 
complications  in  synchronization. 

2.  Spatially  speaking,  for  one-  or  two-dimensional  signal  recordings  within 

a  duration  this  short,  we  should  solve  the  problems  of  light  source  luminosity  and 
rate  of  phenomena  recording. 

3.  The  entire  layout  is  in  an  environment  of  high-tension  discharge  puncturing 
the  air,  ionization  and  interference  of  high-energy  electromagnetic  waves.  For 
example,  the  discharge  voltage  of  a  CO,,  pulsed  laser  is  40-50  kV  and  an  initiating 
pulse  of  40-50  kV  is  the  output  of  high-speed  camera  control  platform.  These 
phenomena  seriously  affect  normal  operation  of  instruments  and  circuits,  causing 
difficulties  in  experiments  and  studies.  51 


4.  A  display  device  should  be  used  to  investigate  non-instantaneous  illumina¬ 
tion,  propagation  of  shock  waves  in  a  transparent  medium,  and  the  motion  of  a 
transparent  vaporized  matter.  The  light  source  of  the  display  device  should  have 
good  interference  properties  and  meet  the  requirements  of  the  highspeed  camera.  In 
other  words,  the  light  source  should  simultaneously  exhibit  interference  and  high 
luminosity  properties  in  order  to  operate  in  a  pulsed  fashion,  thus  preventing 
overlapping  exposure  film. 

5.  The  layout  should  be  optically  matched.  During  joint  operation  of  the 
interference  instrument  and  high-speed  camera,  there  should  be  optical  matching. 
Otherwise,  the  operation  will  fail. 

Solution  of  Technical  Problem  Areas 

1.  Time  synchronization 

With  the  CO.,  pulsed  laser  used  in  the  experiment,  the  output  light  pulse  lasted 
1  microsecond;  the  half  width  of  the  peak  value  is  100  nanoseconds  while  the  laser 
is  charged  with  a  capacitor  and  the  light  emerges  during  the  discharge  from  an 
initiating  sphere  gap.  The  time  delay  from  the  puncture  of  sphere  gap  to  light 
emergence  is  measured  approximately  in  microseconds.  The  lead  time  of  the  initi¬ 
ating  signals  as  output  of  the  controllable  type  high-speed  camera  is  adjusted  by 
the  lead  angle  of  the  revolving  lens.  However,  the  lead  time  is  only  within  one 
revolution  of  the  revolving  lens.  Thus,  the  emergence  of  phenomena  studied  and 
the  position  of  the  revolving  lens  are  synchronized.  However,  under  special 
conditions,  the  phenomena  studied  may  occur  after  the  revolving  lens  gives  an 
initiating  signal  and  arrives  at  a  working  position  of  the  time-delay  rotation. 

If  only  the  joint  operation  of  the  high-speed  camera  and  the  CO^  pulsed  laser  is 
involved,  their  synchronization  is  simple,  requiring  only  initiating  the  sphere 
gap  by  the  camera.  However,  when  another  synchronization  segment  is  required 
and  this  segment  also  occurs  within  a  short  duration,  this  kind  of  synchronization 
is  relatively  complex.  The  light  source  for  the  interference  instrument  is  a 
He-Ne  laser,  which  is  small  as  the  laser  proper,  and  its  power  source  and  the 
laser  have  good  interference  properties.  The  He-Ne  laser  cannot  operate  in 
pulses  [1];  the  laser  has  a  relatively  long  time  delay  from  the  occurrence  of 


the  initiating  signal  to  light  emergence.  The  laser  has  relatively  a  long  time 
deviation,  about  20  milliseconds;  this  is  the  case  when  considerable  excess 
pressure  builds  up  when  the  initiating  signal  is  to  occur  and  its  circuit  is 
easily  interfered  with.  To  get  accurate  synchronization  and  a  controllable  pulsed 
light  source,  if  the  plan  of  igniting  He-Ne  laser  involves  an  initiating  signal, 
no  time  synchronization  can  be  obtained.  We  suggested  a  plan  of  reverse -program 
operation  by  utilizing  the  operational  properties  of  the  He-Ne  laser.  When  the 
applied  current  is  weaker  than  the  working  current,  the  laser  fails  to  glow.  When 
the  applied  current  exceeds  the  working  current  by  a  certain  amount,  the  laser 
ignites  but  no  light  will  emerge.  If  the  laser  is  so  arranged  that  the  excess- 
current  operation  occurs  only  for  a  short  duration,  no  light  then  occurs.  After 
a  synchronized  signal  is  fed,  and  the  excess  current  rapidly  drops  to  the  working- 
current  level  as  controlled  by  an  electron  circuit,  at  that  time,  light  emerges 
instantaneously  and  the  current  cutoff  is  delayed  for  a  certain  duration,  this 
forms  a  single  light  pulse  with  short  time  deviation  and  a  steady  delay  time.  In 
the  circuit,  controllable  silicon  is  used  as  the  He-Ne  laser  switches.  When  the 
controllable  silicon  switches  are  closed,  an  excess  current  is  adjusted  in  opera¬ 
tion,  and  an  external  signal  cuts  off  the  controllable  silicon  switches,  the 
He-Ne  laser  discharges,  relying  on  the  electric  energy  stored  in  the  capacitor 
of  filter  circuit.  The  current  index  in  the  laser  drops.  At  the  working  current, 
the  laser  liberates  light.  The  light  extinguishes  when  the  current  reaches  a 
certain  threshold  value.  Adopting  this  plan  lowers  the  accuracy  of  time  synchroni¬ 
zation  from  more  than  20  milliseconds  down  to  10  microseconds.  Thus,  the  contin¬ 
uous-wave  He-Ne  laser  is  modified  into  a  controllable  single-pulse  type  light 
source  with  high  accuracy  in  time  synchronization;  once  the  pulsed  light  source 
is  matched  with  a  revolving- lens  scanning  high-speed  camera,  the  time  synchroniza¬ 
tion  between  the  two  is  carried  out. 

2.  Luminosity  of  the  light  source 

Light  sources  in  high-speed  photography  involve  a  certain  luminosity  for 

normal  film  exposures.  Visible  light  at  6328A  as  the  output  of  the  He-Ne  laser 

-2  2 

can  produce  an  energy  level  of  10  erg/cm  [7]  on  a  film  sensitive  to  red  light. 

If  the  interference  elements  of  the  (flat-crystal  type  laser)  dislocation  inter¬ 
ferometer  is  a  piece  of  plane  crystal  without  coating  films  at  both  surfaces, 
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the  light  efficiency  is  only  8  percent.  Efficiency  this  low  is  disadvantageous  to 
high-speed  photography.  In  this  case,  the  minimum  exposure  time  is  calculated 
according  to  the  following  equation 


10-’  (sec)  [7] 

In  the  equation,  t  is  the  exposure  time  in  seconds;  S  is  the  sensitivity 
(ergs/cm  )  of  the  film;  A  is  the  illumination  area  in  cm  ;  and  P  is  the  power 
(watts)  of  the  laser  light  source. 

For  example,  an  1-milliwatt  He-Ne  laser  is  used  as  the  light  source  and  a 
plane-crystal  interferometer  without  film  coating  is  the  main  device.  When  the 
illumination  area  is  dlO  millimeters,  the  minimum  exposure  time  is: 

t  =39.3*10  ^  sec 

To  adopt  to  even  higher  exposure  rates  and  to  shorten  the  minimum  exposure 
time,  we  adopted  a  plane-crystal  interference  element  coated  on  both  surfaces  with 
adequate-reflectivity  layers;  light  efficiency  can  be  quintupled  and  a  minimum 
exposure  time  of  7.8*10  ^  second  can  be  realized.  As  proven  by  experience,  wear 
occurs  in  the  present  optical  system  because  of  the  presence  of  many  optical 
elements.  However,  when  Hangwei  [transliteration] -II  type  film  is  used,  a  normal 
density  can  be  obtained  in  the  film  for  an  exposure  time  longer  than  5  microseconds 
If  a  higher-power  device  is  adopted,  no  doubt  a  shorter  minimum  exposure  time  can 
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be  obtained.  For  a  10-milliwatt  device,  a  minimum  exposure  time  of  7.8*10  sec¬ 
ond  can  be  obtained.  If  a  watt  level  argon-ion  laser  is  used,  a  minimum  exposure 
-8  -9 

time  of  10  to  10  second  can  be  obtained. 

3.  Restraint  of  interference 

To  screen  disturbance  from  the  high-energy  electromagnetic  waves,  a  measure 
of  strong  signal  and  high  attenuation  is  adopted,  and  even  the  He-Ne  laser  has  to 
be  screened.  Because  when  an  intensity  electric  disturbance  signal  is  applied 
to  both  terminals  of  an  He-Ne  laser  with  light  emergence,  the  laser  may  be 
quenched,  or  the  light  intensity  abruptly  fades.  In  ways,  interference  is  restrain 
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4.  Optical  matching  of  high-speed  interference  photography 

The  high-speed  camera  and  interferometer  are  two  independent  optical  instru¬ 
ments.  When  they  are  operated  together,  the  optical  matching  problem  must  be 
solved.  The  effective  visual  field  of  the  interferometer  should  have  an  image  of 
specified  size  formed  by  the  camera.  Thus,  the  visual  field  can  be  adequately 
utilized.  This  should  require  an  attached  lens  to  be  set  in  the  light  channel 
to  achieve  optical  matching  between  the  high-speed  camera  and  interferometer.  A 
parallel  light  beam  with  a  <$120  output  diameter  from  the  interferometer  is  easily 
adaptable  to  combined  operation  with  the  single-lens  high-speed  camera  just  by 
removing  the  objective  lens  of  the  camera  and  installing  a  larger-diameter 
(greater  than  the  light-beam  diameter)  objective  lens  with  proper  focal  length. 
Then  the  lens  is  adjusted  to  form  a  clear  image  of  proper  size  on  the  film.  Usual¬ 
ly  it  is  hard  to  find  a  large-diameter  objective  lens,  but  it  can  be  replaced  with 
a  large-diameter  concave  spherical-surface  reflective  mirror  with  good  results. 
Matching  of  the  interferometer  and  a  single-objective-lens  high-speed  camera  can 
produce  photographs  of  the  growing  plume  of  steam  under  laser  irradiation,  and 
high-speed  interference  photographs  of  supersonic  shock  waves  in  a  wind  tunnel; 
these  can  be  examined  to  in  the  reports  [6,  8],  Combined  operation  of  the  inter¬ 
ferometer  and  the  superhigh-speed  camera  with  a  revolving  lens  scanning  frame-by- 
frame  photographs  is  relatively  complex.  Since  the  optical  system  of  a  revolving 
lens  scanning-type  camera  is  an  instrument  with  multi-lens  combination,  its  field 
lens  and  light-alignment  lens  are  strictly  calibrated  and  fixed  during  design 
and  manufacture.  So  matching  between  this  optical  system  and  the  interferometer 
with  large-diameter  parallel  light  as  its  output  demands  a  solution  of  the 
problem  of  light  attachments  in  front  of  the  field  lens  and  light -a ligament  lens. 
Our  conditions  were  as  follows:  the  output  of  a  plane-lens  interferometer  is  a 
parallel  light  beam  with  120  mm  diameter;  after  image  forming,  the  image  on  the 
film  should  have  a  10  mm  diameter;  all  instrument  lenses  should  not  cut  off 
any  portion  of  the  light  beam,  that  is,  the  dimensions  of  light  spots  at  various 
lenses  should  be  smaller  than  the  effective  lens  apertures;  the  light  spot  at 
the  light  gTating  of  the  light-alignment  lens  should  not  be  too  large,  otherwise 
the  exposure  time  will  be  too  long  and  the  overlapping  parts  of  the  film  exposure 
will  be  too  many  between  two  adjacent  frames;  and  the  dimensions  of  light  spots 
at  the  mechanical  shutter  and  revolving  lens  should  be,  respectively,  smaller 
than  the  shutter  diameter  and  the  projection  dimensions  of  the  revolving  lens 
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at  its  working  position.  Only  when  the  above-mentioned  conditions  are  met,  can 
a  complete  visual  field  of  the  interferometer  be  achieved  on  the  film.  In  order 
to  design  a  large-diameter  objective  lens  with  proper  focal  length  and  appropriate 
arrangement  of  the  light  channel,  in  design  a  large-diameter  reflective-type 
objective  lens  should  be  used.  The  objective  lens  is  a  spherical -surface  reflec¬ 
tive  mirror  with  a  980  ami  focal  length,  and  the  focus  falls  at  the  light  axis 

200  mm  in  front  of  the  field  lens  (Fig.  2).  The  known  quantities  are  as  follows: 

the  focal  length  of  the  field  lens  is  110  mm;  the  relative  aperture  diameter  is 
1:4.5;  the  focal  length  of  the  light -alignment  lens  is  56.3  mm;  and  the  width  of 
light  grating  of  the  light-alignment  lens  is  3.6  mm.  In  Fig.  2,  A=53,  B*1S9, 
and  080.  The  calculations  are  as  follows: 

After  focusing  of  the  parallel-light-diameter  concave  lens  M,  the  image  is 
formed  at  200  mm  in  front  of  the  field  lens.  Through  the  field  lens,  an  image  is 

formed  at  a  position  having  a  distance  of  S'  from  the  field  lens 

s'-s£fieH/ts-ffi,ld)-244-2  “ 

In  the  equation,  S  is  200  mm  and  is  the  focal  distance  of  the  field  lens. 


Fig.  2.  Optical  matching  of  inter¬ 
ferometer  and  superhigh-speed  revolving 
scanning  type  camera. 

Key:  (1)  Interference  element;  (2) 
Parallel  light  beam;  (3)  Revolving  lens; 
(4)  Field  lens;  (5)  Light -alignment  lens; 
(6)  Film;  (7)  Superhigh-speed  camera. 


The  position  from  the  light-alignment  lens  is  n*S'-A-B. 
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The  dimensions  of  light  spots  at  various  lenses  and  light  grating  are  as  fol¬ 
lows  : 

(a)  The  dimensions  of  light  spots  at  the  field  lens: 

(b)  The  dimensions  of  light  spots  at  the  reflective  mirror: 

(c)  The  dimensions  of  light  spots  at  the  light-alignment  lens: 

“field" 

Because  of  the  presence  of  light-alignment  lens,  the  distance  between  the 
second  focus  and  the  light-alignment  lens  is  closer  than  in  the  absence  of  the 
light -alignment  lens.  Actually, 

^actual- ^alignment ^  alignment"* 

We  know  from  foregoing  that  the  dimensions  of  image  formed  on  the  film  are 
as  follows: 

^film  ^ ^'^actual^alignment^^actual 

Insert  all  data  into  the  various  equations,  we  obtain:  the  light  spot  at  the 
field  lens,  024.5;  the  light  spot  at  the  lens  surface,  015.9;  the  light  spot  at 
the  light-alignment  lens,  03.21;  and  the  light  spot  at  the  film,  09.4.  All  dimen¬ 
sions  meet  the  requirements.  The  diameter  of  the  concave  reflective-type  objective 
mirror  is  greater  than  120  mm.  After  installation,  the  measured  dimensions  are 
in  agreement  with  the  calculated  values. 

With  ordinary- light  photography,  the  revolving-lens  scanning-type  optical 
system  will  have  very  small  luminous  flux  of  the  instrument  because  of  the  presence 
of  field  lens,  light-alignment  lens,  and  front -mounted  light  grating.  For  normal 
exposures,  the  luminosity  requirements  on  light  source  are  quite  high,  particularly 
so  in  the  superhigh- speed  photography  range.  However,  in  the  present  layout,  the 
field  lens  and  light-alignment  lens  were  coated  with  a  transparency-enhancement 
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film,  and  in  the  light  channel  design  no  part  of  the  light  beam  is  cut  off,  so 
there  is  generally  no  loss.  The  light  efficiency  of  the  instrument  is  better 
and  exposure  time  estimates  are  not  effected  by  the  instrument  types.  This  means 
that  superhigh-speed  photography  can  be  done  with  low-power  (milliwatt  level) 
light  sources. 

Calculations 

After  the  above  technical  problem  areas  were  solved,  the  layout  was  used  to 
obtain  a  series  of  photographs  showing  interaction  between  laser  and  targets. 
Photograph  1  of  the  photograph  section  I  shows  the  concentration  of  a  TEA  C02 
pulsed  laser  beam  on  a  focus,  producing  air  puncturing  and  motion  of  shock  waves. 
Photograph  lb  in  the  same  section  shows  the  plane  waves  reflected  by  spherical 
shock  waves  from  a  spherical -surface  reflective  mirror.  Each  frame  was  exposed 
for  8  microseconds.  The  power  of  the  pulsed  He-Ne  laser  was  1.5  milliwatts. 
Photograph  2  of  the  same  photograph  section  shows  the  puncturing  and  the  jetting 
of  high-energy  vaporized  matter  from  a  solid  target  when  irradiated  with  a 
pulsed  C02  laser.  The  exposure  time  was  5.2  microseconds  for  each  frame;  the 
same  light  source  was  used.  Thus  far,  the  above-mentioned  light  source  has  been 
used  to  take  high-speed  interference  photographs  of  1.3  microseconds  of  exposure 
for  each  frame.  These  recordings  no  doubt  can  provide  qualitative  and  quantitative 
evidence  of  the  interaction  of  lasers  and  matter. 

In  addition,  the  experimental  results  proved  that  this  investigation  technique 
can  be  used  to  show  more  displays  of  high-ppeed  flow  fields,  like  laser  wind  tun¬ 
nel,  explosion,  combustion,  and  plasma.  Photograph  3  of  the  photography  section  I 
shows  an  investigation,  using  this  installation  of  a  light-cavity  flow  field  of 
an  aerodynamic  laser.  We  can  determine  the  effect  on  quality  of  output  light  beam 
by  the  position  of  shock  waves  produced  at  the  throat  channel .  Photograph  4  of  the 
photography  section  I  shows  the  flow-field  photograph  near  a  model  placed  in  a 
M  2.5  wind  tunnel.  The  two  above-mentioned  photographs  were  taken  by  using  an 
1-milliwatt  He-Ne  laser  with  continuous  wave  output  as  the  light  source;  each 
frame  exposure  is  0.5  and  1  millisecond,  respectively. 
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In  this  paper,  we  set  up  a  three-dimensional  statistical 
simulation  model  for  the  problem  of  rarefied  gas  hypersonic  flow 
past  combined  bodies  of  revolution  with  an  angle  of  attack.  In 
order  to  improve  the  efficiency  of  the  simulation  and  to  expand 
the  simulation  region,  we  adopt  a  special  weighting  technique 
and  reduced  collision  model.  A  series  of  calculation  has  been 
executed  on  013  Computer  and  the  results  are  very  satisfactory. 


I .  Foreword 


Studies  of  aerodynamic  properties  in  rarefied  gases  are  decidedly  important 
for  space  flight  vehicles  (extended  flights),  such  as  meteorological  satellites, 
manned  and  other  spacecraft  [1].  In  recent  decades,  researchers  studying  rarefied 
gas  aerodynamics  have  persisted  in  trying  to  solve  the  Boltzmann  equation,  yet  up 
to  now  no  important  breakthrough  was  forthcoming.  For  urgent  and  practical 
problems,  they  rely  only  on  experiments  and  engineering  methods.  However,  the 
statistical  simulation  method  developed  by  Bird  can  provide  almost  all  aerodynamic 
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data  [2,  3]  for  flight  vehicles  over  the  entire  transitional  region.  Relying  on 
Bird's  contribution  [2]  and  the  report  [4],  the  authors  established  a  statistical 
model  of  combined  bodies  of  revolution  for  three-dimensioanl  circulation  flow;  a 
series  of  calculations  were  made  for  a  short  conical  frustum  at  its  base  of  the 
spherical  platform,  to  derive  the  pressure  coefficient  c^,  friction  coefficient 
c^,  heat  transfer  coefficient  st,  coefficient  of  axial  force  c^,  coefficient  of 
normal  force  c^,  moment  coefficient  cn>  pressure  center  position  x  ,  and  other 
parameters.  For  better  simulation  efficiency  and  longer  computation  range,  a 
method  of  simulated  molecular  weighting  is  introduced  based  on  Bird's  approach 
[2].  The  computational  results  are  in  general  agreement  with  the  experimental 
results.  This  paper  briefly  introduces  some  computational  results. 


II.  Description  of  Method 


A  statistical  simulation  method  is  actually  a  gas  flow  process  simulated  by 
a  computer;  a  space  region  is  marked  out  and  referred  to  as  a  flow  field.  This 
paper  adopts  a  cylinder- shaped  flow  field;  an  object  is  placed  in  the  central 
portion  of  the  flow  field;  the  object's  symmetrical  axis  is  coincident  with  that 
of  the  flow  field.  The  incoming  flow  is  placed  in  the  XY  plane.  From  these 
assumptions,  it  is  sufficient  to  only  simulate  one  half  of  the  space  (see  Fig.  1). 
The  study  [4]  describes  the  grid  division,  simulation  of  the  initial  position  of 
molecules,  determination  of  velocity,  and  methods  of  supplementing  molecules  at 
the  inlet  of  the  flow  field.  The  simplification  and  treatment  methods  of  colli¬ 
sions  are  borrowed  from  the  report  [2].  The  molecules  exhibit  a  mirror -surface 
reflection  to  the  symmetrical  plane  (Z=0 ) .  The  reflection  model  on  the  object 
surface  is  cited  from  the  studies  [2,  3]. 

To  extend  the  simulation  range  and  for  better  simulation  efficiency,  this 
paper  adopts  a  weighting  technique,  by  assigning  different  weights  to  simulation 
molecules  of  different  parts  of  the  flow  field.  In  other  words,  weights  represent 
actual  gas  molecules  in  different  numbers.  The  method  used  in  the  paper  consists 
of  assigning  weights  along  the  diametral  direction  of  the  flow  field.  Let  r.  ,  . 

*  f  J  9  * 

stand  coordinates  at  the  diametral  direction  of  grid  nodal  points  with  subscripts 
(i,  j,  k).  Because  of  different  weights,  when  a  simulated  molecule  enters  another 
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grid  from  a  grid,  the  molecule  may  disappear,  continuously  exist,  or  divide  into 
several  simulated  molecules.  For  an  example,  when  a  simulated  molecule  enters 
another  grid  from  a  grid,  the  molecule  may  disappear,  continuously  exist,  or  divide 
into  several  simulated  molecules.  For  an  example,  when  a  simulated  molecule 
enters  another  grid  with  subscripts  (i,  j,  k+1)  from  a  grid  (i,  j,  k),  should  we 
use  a  random  number  to  determine  whether  or  not  the  molecule  continuously  exists? 


When 


Rt< 


±J- 


I'i.i.t+t  +  * 

the  molecule  exists  continuously.  Otherwise,  the  molecule  disappears.  In  the 
equation,  is  a  random  number  with  a  uniform  distribution  in  (0-1)  regions.  If 
a  simulated  molecule  enters  another  grid  with  subscripts  (i,  j,  k+1)  from  a  grid 
(i,  j.  k) ,  then  a  simulated  molecule  may  divide  into  several  molecules;  the 
integer-number  part  of  the  expression  equation  (2r.  .  .  +l)/(2r.  .  .+1)  shows  the 

J  » K+i  1 1 J t K 

number  of  molecules  after  dividing.  If  the  part  following  the  decimal  point  is 
greater  or  equal  to  the  random  number  R.^,  to  the  number  of  daughter  molecules  is 
added  1;  otherwise,  no  addition  is  made.  For  the  situation  of  several  daughter 
molecules,  it  is  necessary  to  artificially  handle  the  molecule  position  and  velocity 
in  order  to  have  a  slight  difference  between  the  velocities  of  daughter  molecules 
and  to  ensure  that  any  point  in  space  can  be  occupied  by  one  simulated  molecule 
at  most.  Thus,  compute  the  collisions  between  molecules  and  trace  the  variation 
in  position  and  velocity  of  each  simulated  molecule.  The  recordings  and  their 
proper  treatment  can  lead  to  various  parameters  representing  aerodynamic  properties. 


III.  Analysis  of  Results 


By  using  the  above-described  method,  in  the  paper  a  series  of  computations 
of  a  combined  body  of  revolution  having  circulation  flow  with  angle  of  attack  is 
made.  The  exterior  of  the  combined  body  of  revolution  is  shown  in  Fig.  2.  The 
variation  range  of  incoming-flow  parameters  is  16-22  (s  is  the  ratio  of  the 
incoming-flow  velocity  and  the  thermodynamic  velocity  of  the  molecules);  the 
ratio  of  wall  temperature  and  total  temperature  tw=0.1  to  0.13;  the  angle  of 
attack  0*0  to  20°.  Computations  were  conducted  on  a  Chinese-made  013  Computer; 
the  operating  speed  of  the  floating  point  mode  of  the  computer  is  2  million  times 
per  second,  and  its  internal  memory  can  hold  130,000  words.  The  number  of  initial 
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molecules,  without  weighting,  are  3*10  to  10  .  The  dimensionless  length  of  time 

\  % 

step  is  At  at  0.01-0.05  ( At= At Adhere  At,  cm00  and  A^  represent,  respectively, 

the  length  of  time  segment,  thermodynamic  velocity  of  molecules  in  incoming  flow, 

and  the  average  collision  freedom).  In  order  to  eliminate  as  much  as  possible  the 

influence  of  non-equilibrium  effect,  sampling  begins  after  a  40-60At  time  interval 

the  sampling  period  is  40-80 At.  At  zero  angle  of  attack,  there  are  two  regions 

along  the  direction  of  the  argument;  there  are  divided  into  three  or  four  regions 

when  the  angle  of  attack  is  not  equal  to  zero.  When  the  weighting  technique  is 

3  3 

used,  the  number  of  initial  molecules  is  l.S-10  -5‘10  . 


Fig.  1.  Division  of  three-dimensional  flow  field  and  grid. 
Key:  (*)  Model. 


Figure  2  shows  a  velocity  diagram  without  weighting;  data  are  obtained  from 
the  first  and  fourth  region  in  the  direction  of  the  argument;  these  can  be  approx¬ 
imately  considered  as  the  windward  and  leeward  surfaces;  the  length  of  line 
segments  stands  for  the  magnitude  of  velocities;  directions  are  represented  by 
arrows.  Since  the  velocity  component  at  Z  direction  is  neglected,  the  diagram 
gives  -+-  v*  as  the  actual  figure.  The  velocity  diagram  shows  some 
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Fig.  2.  Velocity  spectrum  of  (V^+Vy)1^2 

=-0.167,  6-20 

*D«  ' 

fw  *»  0.1,  a  «*  10* 


flow  details.  First,  near  the  stationary  point,  the  magnitude  and  direction 
of  velocity  have  changed  apparently.  Secondly,  the  disturbance  at  the  windward 
surface  of  the  object  is  greater  than  the  leeward  surface.  Lastly,  from  an  analysis 
of  the  flow  field,  we  can  see  that  when  the  angle  of  attack  a  is  increased,  the 
diametral -direction  length  of  the  flow  field  also  correspondingly  increases  in 
order  to  compare  two  reflection  models  from  the  wall  surface  and  to  calculate 
the  situation  of  two  sets  of  zero  angle  of  attack.  Fig.  3  gives  the  comparison 
between  the  calculated  drag  coefficients  and  experimental  results  [5].  Basically, 
the  calculated  data  are  in  agreement  with  the  experimental  results.  There  are 
little  differences  in  effects  between  two  models;  this  is  possibly  because  the 
incoming  velocity  is  quite  high. 

Variations  of  aerodynamic  coefficients  Cy  c^,  and  cm  with  the  angle  of  attack 
a  are  given  in  Fig.  4.  a  is  0°,  5°,  10°,  and  20°.  In  these  figures,  wheno=0°, 
the  sampling  period  is  4Q£tt;  in  other  cases,  the  sampling  period  is  8C&^.  The 
darkened  symbol  when  ct=108  represents  the  computational  result  without  weighting. 

It  is  apparent  that  the  weighting  method  is  successful.  From  the  figure,  we  also 
can  see  linear  variations  of  c^  and  cm  with  the  angle  of  attack;  the  result  is 
completely  identical  with  the  conclusion  obtained  from  the  rarefied-gas  experiment 
when  the  angle  of  attack  is  small. 
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A  Fix:  Tunnel  *L» 
OE  Shock  Tunnel 
GE  Shock  Tunnel 
C if  Shock  Vurtnel 
VKL  Fret  Jei 
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Fig.  3.  Comparison  between  theoretical  values  and 
experimental  values  (drag  coefficient  c_). 

Key:  (a)  Symbols;  (b)  Layout;  (c)  Semi-conical  angle; 
(d)  This  paper;  (e)  Reflection  model  in  12];  (f)  Reflec 
tion  model  in  [3];  (g)  (Obtuse). 


i 
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This  paper  derives  an  isothermal -curve  diagram  (Fig.  5)  concerning  tempera¬ 
tures;  the  curves  connect  various  points  of  same  temperature  in  the  flow  field. 

Since  the  wall  temperature  is  relatively  low,  at  the  nose  portion  a  high-temperature 
region  is  formed;  the  curve  is  generally  symmetrically  distributed  near  the 
stationary  point. 


0  Cj 

0  CM 

&  C. 


Fig.  4.  Variations  of  c^,  c^, 

and  c  with  a. 
m 

O  CA  □  cH  A  cm 


IV.  Conclusions 

In  the  paper,  a  series  of  computations  are  conducted  concerning  rarefied-gas 
circulation  flow.  As  indicated  by  results,  a  three-dimensional  statistical  model 
is  established  on  the  basis  of  two-dimensional  and  symmetrical  circulation  flow 
calculations  as  reported  by  the  studies  [2-4];  both  the  three-dimensional  statis¬ 
tical  model  and  the  adopted  weighting  method  are  successful. 
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